Numerical study on shock–dusty gas cylinder interaction

Jingyue Yin·Juchun Ding·Xisheng Luo·Xin Yu

Abstract The interaction of a planar shock wave with a dusty-gas cylinder is numerically studied by a compressible multi-component solver with an adaptive mesh refinement technique.The influence of non-equilibrium effect caused by the particle relaxation,which is closely related to the particle radius and shock strength,on the evolution of particle cylinder is emphasized.For a very small particle radius,the particle cloud behaves like an equilibrium gas cylinder with the same physical properties as those of the gas–particle mixture.Specifically,the transmitted shock converges continually within the cylinder and the n focuses at a region near the downstream interface,producing a local high-pressure zone followed by a particle jet.Also,noticeable secondary instabilities emerge along the cylinder edge and the evident particle roll-up causes relatively large width and height of the shocked cylinder.As the particle radius increases,the flow features approach those of a frozen flow of pure air,e.g.,the transmitted shock propagates more quickly with a weaker strength and a smaller curvature,resulting in an increasingly weakened shock focusing and particle jet.Also,particles would escape from the vortex core formed at late stages due to the larger inertia,inducing a greater particle dispersion.It is found that a large particle radius as well as a strong incident shock can facilitate such particle escape.The theory of Luoetal.(J.Fluid Mech.,2007)combined with the Samtaney–Zabusky(SZ)circulation model(J.Fluid Mech.,1994)can reasonably explain the high dependence of particle escape on the particle radius and shock strength.

Keywords Dusty-gas cylinder·Non-equilibrium effect·Shock wave·Instability

1 Introduction

Dusty-gas flows,which refer to flows of gases seeded with small solid particles,occur widely in various industrial applications such as explosive safety[1] and plasmas of the solar system[2].In dusty-gas flows,particles and gas exchange their momentum and energy continuously in the transitional zone, and this non-equilibrium effect would significantly modify the fluid properties and further the flow characteristics.For instance,in advanced flow diagnostic techniques such as particle image velocity(PIV)[3] and planar laser Rayleigh scattering[4],a large number of small particles are injected into the test gas serving as tracers.If the radius or mass fraction of the particles is very small,the resulting image quality will be poor,whereas if it is relatively large,the solid phase would have a considerable interaction with the gas flow,causing a large deviation of the observed result from the actual flow.

Shock–cylinder interaction is a fundamental configuration in studying shock-accelerated interfacial instability[5–7].When a shock wave passes over a cylindrical interface separating two materials of different properties(e.g.,density,specific heat ratio),complex flow phenomena emerge including the shock wave refraction,the generation and transport of baroclinic vorticity, and continual deformation of the cylindrical interface.Numerous the oretical models on shock–cylinder interaction have been proposed, and some of them have been validated against experiments and simulations.The deformation of a spherical or cylindrical gas inhomogeneity was first proposed by Rudinger and Somers[8] and a theoretical model was established to predict the interface velocity.Later,interactions of a planar shock wave with heavy and light gas cylinders or bubbles were experimentally investigated, and the wave configurations as well as the velocities of the shock waves and interfaces were analyzed[9].Also,advanced techniques such as the planar laser induced fluorescence,the Mie scattering, and the high speed Schlieren photography were applied to capture the detailed process of shock–gascylinder interaction under various Mach numbers[10,11],gas species[12,13] and aspect ratios[14].Recently,the evolution of three-dimensional light and heavy gas cylinders with controllable initial shapes impacted by a weak shock wave was studied both experimentally and numerically[15,16].Most of previous studies focused on the development of cylindrical interface between two materials of the same phase,but this is not the case for real applications such as explosions with reactive metal particles[17]where the materials on both sides of the interface usually possess different phases.

Few works on shock–dusty gas interaction have been reported[18,19].The deformation of a dense-gas region subjected to a planar shock was studied by Otaet al.[20], and a“roll-up”structure together with a“mushroom”jet was definitely formed.A similar setup using particle clouds instead of heavy gas was simulated by Kiselev et al.[19], and the particle evolution was found to be similar to that of gas cylinder.The particle radius in the se works is very small, and ;thus,the non-equilibrium effect associated with the particle relaxation was negligible.Later,the interaction between a planar shock and a dusty-gas single-mode interface was investigated by a Lagrangian tracking method, and a linear model for two-phase Richtmyer–Meshkov(RM)instability was proposed[21].It was found that for duty gases with the Stokes number less than 1.0,which is heavily dependent on the particle radius,mass fraction, and shock strength,the growth of the perturbation amplitude was very close to the equilibrium gas case,whereas if the Stokes number is much larger than 1.0,the linear growth mode would convert to be exponential.Also,the interaction between a planar shock and the dusty-gas was experimentally studied by Vorobieff et al.[22].The physical mechanism for the perturbation growth of two-phase RM instability was found to be quite different from that of the single-phase case due to the existence of non-equilibrium effect caused by particle relaxation.

The non-equilibrium effect behaves differently for dustygas flows with different characteristic times.Specifically,when the timescale of particle relaxation is much larger than the flow characteristic time(e.g.,for the case of a large particle radius and small flow velocity),the flow is close to a frozen flow,whereas if the relaxation time is comparable to the flow characteristic time,the non-equilibrium effect plays an important role on the flow structure.For the latter,taking the shock propagation in dusty gases as an example,across the leading front of a dusty shock,the gaseous phase experiences a sudden change in its flow property.For the larger inertia,the particles could not follow the rapid change in gas velocity and temperature,exhibiting a relaxation region behind the shock front.Through out the transition zone,particles and gas exchange their momentum and energy continuously, and finally an equilibrium state was established at the end of the relaxation zone[23].For the shock reflection over a single wedge in a dusty gas,an abnormal transition between regular reflection and Mach reflection was observed, and the wave structure propagating along the wedge surface was also found to be non-self-similar[24].For the dusty shock reflection over a double wedge,the length of the former wedge influences the reflection the back wedge, and some reflection types different from those in the pure gas case were observed[25].It is seen that little attention has been paid to the non-equilibrium effect on shock–dusty gas cylinder interaction,which motivates the present work.

In this paper,the evolution of a cylindrical particle cloud impacted by a planar shock is studied using a compressible multi-component solver.The dusty-gases with different particle radius(1,2,4, and 8μm)impacted by a planar shock with different strengths(Ms=1.2,1.65, and 2.0)are considered.We mainly discuss the initial condition dependence on the flow structures including the density distributions of gas and particles,the wave patterns,the pressure contour, and the interface distortion.Also,a detailed comparison for the perturbation growth rate between the numerical results and the two-phase model of Ukai et al.[21]is given.

2 Numerical and theoretical overview

A series of assumptions are adopted in the present simulations,which are also commonly accepted in previous works:

(1)solid particles are uniformly distributed in the carrying gas;

(2)the particle volume as well as their collision is not considered;

(3)the number density of particles(number of particles per unit volume)is large enough so that the particles could be considered as a continuum;

(4)particles are spheres with the same attributes including the diameter,density,specific heat ratio;

(5)the gas is treated as ideal and the viscosity and heat transfer within the gas phase are ignored;

(6)the transport process between the gas and particles is modeled by drag forces and heat conduction.

Dusty-gas flows in two-dimensional geometry can be described by the conservation laws of mass,momentum and energy for each phase as follows:

where

Here ρ,p,E,u, and v denote the density,pressure,specific total energy, and velocity components in the x and y directions,respectively.The subscripts g and p indicate the gaseous and solid phases,respectively.The perfect equation of state is employed and the energies for gas and particles are respectively expressed as:

where T,Cv, and Cmst and for the temperature,the specific heat capacity of gas at constant volume, and the specific heat capacity of particle material,respectively.In the current study,the detailed transport process between particles and gas is treated independent of the flow dynamics and characterized by the drag forces(Dx,Dy) and heat transfer(Q),as included in the source terms of the control equations.The y are described as:

where ρdis the material density of particle,rdis the particle radius,Cpthe specific heat capacity of gas at constant pressure,Nuthe Nusselt number, and Prthe Prandtl number.

It was found that for low mass loading ratios,the influence of early-stage variation in drag force was negligibly small[26], and thus the steady drag coefficient model is employed in the present simulations.Several steady models[27–29]were examined, and the calculated shock profiles agree well with each other,demonstrating a negligible influence of the drag coefficient on the computational results.In this work,the coefficient model of Crowe[27]is used for all simulations.The dependence of the Nusselt number model on the shock wave profile was checked by Saito et al.[24], and the calculated shock profile employing the Knudsen and Katz model[30]agrees well with other models[31–33].In this work,the simple but effective Knudsen and Katz model is employed,which is expressed as Nu=2.0+0.6P1/3rR1/2e,with Rebeing the particle Reynolds number.

Fig.1 a Evolution of the shocked single-mode dusty interface and b the corresponding amplitude growth under different initial particle mass loadings with a particle radius of 0.5μm.The black lines indicate the results of Ukai et al.[21], and the symbols represent the present numerical results

The governing equations are solved by a finite volume method based on the unstructured mesh discretization.Numerical fluxes at the cell edges are calculated with the MUSCL-Hancock scheme, and second-order accuracy is achieved in both time and space.The present solver has been widely used in two-phase flow simulations[25], and showed good performance.To further examine the solver capacity in resolving dusty-gas flows,the problem of a shock wave interacting with a single-mode dusty-gas interface by Ukai et al.[21]is simulated.The corresponding computational configuration and the subsequent interface evolution are shown in Fig.1a.After the incident shock impacts the single-mode interface,a reflected and a transmitted shock are immedi-ately produced.Also,the interface deforms gradually with an increasing amplitude, and later nonlinear bubble(light fluid penetrating into heavy one) and spike(heavy fluid penetrating into light one)structures are formed.Figure 1b presents the amplitude growth of the shocked interface for four typical cases with different particle mass fractions but the same particle radius(0.5μm).As we can see,good agreement is obtained between the present simulations and previous results for all cases,which demonstrates the good reliability of the present solver.

The numerical configuration corresponding to the interaction of a planar shock with a dusty-gas cylinder is shown in Fig.2,where the planar shock propagates from left to right with various strengths of Ms=1.2,1.65, and 2.0.The length and height of the computational domain are 200 mm and 70 mm,respectively.The left and right edges are treated as flow-through boundary so that the outgoing waves can exit the domain without reflection, and the other boundaries are set as symmetric or reflecting.The direction along and perpendicular to the symmetric boundary are set as x and y axes.The dusty-gas cylinder with an initial radius of R=17.5mm is located 15mm downstream from the shock wave.The solid particles with radii of 1,2,4, and 8μm are uniformly distributed in the gas cylinder surrounded by the environmental air.For different cases,the mass loading is fixed and a constant Atwood number of At=0.639 is maintained.Detailed properties of the gas and particles are listed in Table1.The pressure and temperature of the environmental air are pg=1 at m and Tg=298 K,respectively.The material density of particles is ρd=780 kg/m3and the particle temperature is Tp=298 K.The particle phase density is ρp=3.535 ρg,the specific heat capacity of air at constant volume is Cv=0.717 kJ/(kg·K), and the specific heat capacity of particles is Cm=1.726 Cv.For particles with a very small radius(1 μm),the gas–particle mixture inside the cylinder be haves like an equilibrium heavy gas with effective physical properties.To highlight the particle relaxation effect on the inter facial evolution,the dusty-cylinder evolution is compared with that of equilibrium gas case under the same initial conditions.

Fig.2 Schematic of the initial computational setting

Table1 Initial gas properties inside and outside the cylinder,where the interior dusty gas is treated as a heavy equilibrium gas

Fig.3 Density distribution along the x-axis at 200μs under different mesh sizes,where the pre-shock air density(ρ0)is chosen as a characteristic value

A mesh refinement study is performed to ensure the simulation accuracy.The density profiles for a gas–particle mixture with a particle radius of 1μm under different mesh sizes of 1,0.5,0.25, and 0.125mm along the x-axis are given in Fig.3.It is shown that the numerical solution converges as the mesh size reduces from 1 to 0.125 mm.Since the present work mainly focuses on the intermediate-to large scale structures of the interface,a compromise between the solution accuracy and the computational cost can be made.Hence,the initial mesh size of 0.25 mm is employed in the present simulations.To effectively capture the wave patterns and density inhomogeneities,an adaptive mesh refinement technique is employed.The adaptation procedure refines the quadrilateral mesh in flow regions with large truncation errors and coarsens the grid in regions with small truncation errors,while the basic mesh is always retained.A measure for the truncation error at a cell edge is given by the ratio of the second derivative and the first derivative in the Taylor series expansion of the density.If this truncation error is larger than a specified threshold for refinement,the cell will be refined.While if this error is smaller than a specified threshold for coarsening,the mesh will be coarsened.To give the readers a clear image,a refined mesh for the shock–dusty–gas cylinder interaction at 1400μs is displayed in Fig.4.As we can see,the present adaptive mesh refinement technique can well refine the flow regions with large density gradient.Benefiting from this,we are able to obtain numerical results of considerable accuracy under acceptable computational sources and time.The zero time is defined as the moment when the incident shock meets the leftmost boundary of the cylinder.

Fig.4 Refined mesh for the shock–dusty–gas cylinder interaction at 1400μs

3 Results and discussion

The dependence of cylinder evolution on the particle radius,which is a key parameter influencing the non-equilibrium effect,under a Ms=1.2 shock wave is first considered.The density distributions of particle and gas,wave patterns, and pressure contours at the early stage are discussed in detail.Temporal variations of the leftmost boundary position and the dimensions of the evolving cylinder(width and height)are also given.In addition,special attention is paid to the effect of shock strength on the behavior of particle escape from the vortex core.

3.1 Density contour and pressure profile

The density distributions of gas(color flood) and particles(black line)for the interaction between an Ms=1.2 planar shock and a dusty-gas cylinder with varying particle radius(rd=1,2,4, and 8μm)are shown in Fig.5.Complex wave structures are formed at the early stage of shock–cylinder interaction.As the incident shock hits the upstream interface of the cylinder,a curved transmitted shock propagating in dusty air and a reflected shock moving upstream are formed.For particles with a small radius(1μm),the transmitted shock converges continually and the n focuses at a small region within the cylinder,producing a local high-pressure zone(160μs).Since this high-pressure zone is very close to downstream interface,a downstream-moving particle jet is subsequently induced(400μs).Here,the fastest-developing jet occurs at a position slightly above the symmetry axis.This is different from the results of shock–heavy gas cylinder interaction where the fastest jet usually develops along the symmetry axis.A possible reason is that particles along the symmetry axis have a larger number density caused by the stronger focusing there and ,hence,attain a lower speed under a constant driving force.As the particle radius increases to 2μm,the jet develops much more slowly due to the weaker shock focusing.For particle radius lager than 2μm,generally,the flow features approach those of a frozen gas flow of pure air,e.g.,the transmitted shock increases in velocity but decreases in curvature.Particularly,for the largest radius case,the gas-dust cylinder only suffers the shock compression without any distortion.

The pressure distributions along the x-axis at two typical instants are plotted in Fig.6,where the results of equilibrium gas case are also given for comparison.The color lines with symbols indicate the dusty gas with different particle radius, and the black line represents the equilibrium gas case with the same physical properties as those of gas–particle mixture listed in Table1.The pressure distribution for the 1μm case at 20μs is similar to the equilibrium gas case, and there exists a rapid pressure jump across the transmitted shock front.As the particle radius increases,the transmitted shock propagates more quickly with a weaker strength and a smaller curvature.More importantly,the pressure rises after the transmitted shock front becomes more gentle,which indicates a longer relaxation zone.Particularly,for the 8μm case,the pressure suffers a sudden jump across the transmitted shock front, and the n increases gradually through the relaxation zone where particles and gas exchange their momentum and energy and finally reach a balance.As the shock focusing occurs,a pressure peak is produced as shown in Fig.6b.As we can see,the pressure peak for the 1μm case is nearly the same with the equilibrium case due to the weak non-equilibrium effect.When the particle radius is larger than 2μm,the pressure peak drops significantly for the much weaker shock focusing.This provides a strong evidence for the disappearance of jet structure in large particle radius cases.It is interesting that for all dusty cylinders,there exits a second pressure peak in front of the transmitted shock.This is produced by the reflection of the diffracted shock from the x-axis where a reflecting boundary condition is imposed.As the particle radius increases,the diffracted shock becomes weaker, and ;thus,the second pressure peak drops and even disappears in the 8μm case.Note the transmitted shock in the equilibrium gas case is much slower than dusty-gas cases, and ;thus,no second pressure peak emerges at the early stage when the diffracted shock has not reached the x-axis.

3.2 Interface deformation

Schlieren results corresponding to the density gradient of gas–particle mixture for cases of different particle radius are shown in Fig.7.Since the density gradient for gas phase is very small,the present Schlieren results mainly reflect the density variation of particles.Numerous interesting flow phenomena are observed.The particle jet is formed at early stages for small particle radius cases,which is similar to the equilibrium gas case.The main reason is that particles with a small size have a small inertia and thus can follow the gas flow quickly,i.e,a very weak non-equilibrium effect.In addition,for the small radius cases,the re exist noticeable secondary instabilities along the shocked cylinder edge at late stages, and the particles are also able to“roll up”due to the small inertia.For the larger radius cases,particles would escape from the vortex core formed at late stages,causing a greater dispersion of particles(i.e.,larger width and height of the particle cylinder).The reason is that the negative pressure gradient with in the vortex core cannot support the centripetal force of the large-inertia rotating particles.This is more evidently shown in Fig.8 where the density contours of gas and particles with rd=1 and 8μm at 700μs are compared.For the 1μm case,the high gas density region also has a high particle density,which indicates the good following perfor-mance of particles.For the large radius case,the particles aggregate at the upper right of the vortex core.

Fig.5 Density contours for gas and particles at typical moments.IS:the incident shock,TS:the transmitted shock,DS:the diffracted shock

Fig.6 Pressure profiles along the x-axis for dusty gas with different particle radius at a 20μs and b 60μs

Fig.7 Typical Schlieren pictures corresponding to the density gradient of gas–particle mixture for a M s=1.2 shock interacting with dusty-gas cylinders with different particle radius

Fig.8 Gaseous and particle density contours with particle radius a 1μm and b 8μm at 700μs

Previous work of Ukai et al.[21]found that a shocked dusty-gas perturbed inter face suffers a linear growth in amplitude for Stokes number less than 1.0(similar to that of equilibrium gas case),but increases exponentially for Stokes number larger than 1.0.The Stokes number is defined as

Temporal variations of the height and width of the semis hocked cylinder are given in Fig.10.For cases with a particle radius smaller than 2μm,the cylinder dimension variation is similar to that of heavy gas cylinder[16].It is found that the width of the shocked cylinder for small radius cases is larger than that of large radius cases at the early stage.The main reason is that the “roll-up”near the uppermost point of the cylinder in smaller radius cases greatly distorts the interface.Because particles with larger inertia are easy to depart from the gas flow,the shocked cylinder filled with large-sized particles tends to have larger width and height.

Fig.9 Trajectory of the left most boundary under different particle radii

Table2 Velocity of the leftmost boundary of a cylinder as well as the Stokes number.The prediction of Wang’s theoretical model[34]is given for comparison

3.3 Particle escape from vortex core

In this section,we mainly study the dependence of particle escape behavior on the shock Mach number.The interaction between a planar shock wave of different strengths(Ms=1.65 and 2.0) and a dusty-gas cylinder with various particle radii(rd=1,2,4,8μm)is considered.

Representative Schlieren images illustrating the detailed deformation of gas–particle mixture for all cases are displayed in Fig.11.For stronger incident shocks,the particle jet as well as the secondary instability emerges only if the particle radius is very small,which is similar to the weak shock situation shown in Fig.7.A comparison among cases of different shock strengths shows that the larger the Mach number,the stronger the secondary instability, and thus the greater the cylinder deformation.It indicates that,in addition to the large particle radius,a strong shock wave can also facilitate the particle escape from the vortex core.Taking the 2 μm case as an example,particles can well follow the rotating air near the vortex core under Ms=1.2,but largely escape from the vortex core for Ms=1.65 and 2.0.It means that a stronger shock can induce a greater dispersion of particles, and thus larger width and height of the shocked cylinder.

Fig.10 Time-variation of a the height and b the width of the shocked cylinder under different particle radius

Fig.11 Density gradient of gas–particle mixture for a planar shock of M s equals a 1.65 and b 2.0 interacting with dusty-gas cylinders with different particle radius

In previous work of Luo et al.[36],the onset time of particle escape from a vortex core τdwas found to rely on the particle radius rdand the vortex circulation Γ, and can be estimated as where Ckdenotes the Knudsen correction for the Stokes law, and usually takes a value of 0.6.The ω refers to the local vorticity, and rcto the vortex radius,which can be measured by the distance of the vortex center from the uppermost point of the shocked cylinder.

Previous experiments on shock–heavy gas cylinder interaction have observed a pair of counter-rotating vortices at the late stage due to the evolution of baroclinic vorticity initially deposited on the density interface[5].Several theoretical models were proposed to calculate the vortex strength[8,37,38],among which the Samtaney–Zabusky(SZ)model is considered to be the most accurate one due to its strict theoretical derivation based on shock polar analysis.It can be expressed as

Table3 Particle escape time predicted by Eqs.(3) and (4)

where a is the sound speed.

Since the circulation deposition is mainly determined by the pressure gradient of the incident shock and the density gradient across the material interface,the SZ model can also give a reasonable estimation of the circulation deposited during the present shock–dusty gas cylinder interaction.Both the model prediction and the present simulation show that the vortex circulation Γ augments with increasing the shock strength.The combination of Eqs.(3) and (4)can give a quantitative prediction of τdwith the known particle radius and shock strength.As shown in Table 3,the onset time of particle escape for the present shock–dusty gas cylinder interaction decreases continually with increasing either the particle radius or the shock strength.It illustrates that particles are easier to escape from the vortex core under greater shock strength and larger particle size.The theoretical prediction reasonably explains the high dependence of escape time on the particle radius and shock strength found by simulations.Note the present findings are only applicable to shock–dusty gas cylinder interactions under small particle radius,low mass loading ratio, and weak shock strength.Further study considering a wider range of parameter spaces will be reported in the near future.

4 Conclusion

In this work,the interaction between a planar shock wave of different strengths and a dusty-gas cylinder with varying particle radii has been studied by an adaptive compressible multi-phase solver.Special attention was paid to the particle radius influence on flow structures including the wave propagation,the jet formation, and the cylinder deformation.For a very small particle size,the particle cloud behaves like an equilibrium gas cylinder with the same physical properties as those of the gas–particle mixture.Specifically,the transmitted shock converges continually within the cylinder and the n focuses at a region near the downstream interface,producing a local high-pressure zone,followed by a particle jet.The re exist noticeable secondary instabilities along the cylinder edge and the evident particle roll-up causes relatively large width and height of the shocked cylinder at the early stage.As the particle radius increases,the transmitted shock propagates more quickly with a weaker strength and a smaller curvature, and thus the shock focusing and the particle jet become much weaker and even disappear.It demonstrates that with increasing the particle radius,the flow features approach those of a frozen gas flow of pure air.Moreover,particles would escape from the vortex core formed at late stages due to the large inertia,causing a greater dispersion of particles(i.e.,larger width and height of the particle cylinder at late stages).

The effect of shock strength on the particle escape behavior is also discussed.It is found that the larger the Mach number,the stronger the secondary instability, and thus the greater the cylinder deformation.In addition to the large particle radius,the strong shock wave can also facilitate the particle escape from the vortex core.The high dependence of particle escape on particle radius and shock strength can be reasonably explained by the theory of Luo et al.[36]combined with the SZ circulation model.The present findings are useful for assessing the effect of tracing particles on the actual gas flow in advanced flow diagnostic techniques such as PIV.

AcknowledgementsThis work was supported by the National Natural Science Foundation of China(Grants11802304 and 11625211) and the Science Challenging Project(Grant TZ2016001).