1689881

ContentslistsavailableatScienceDirect

International Journal of Multiphase Flow

journalhomepage:www.elsevier.com/locate/ijmulflow

Computational analysis of shock-induced flow through stationary particle clouds

Andreas Nygård Osnes

^a,∗

, Magnus Vartdal

^b

, Marianne Gjestvold Omang

^c,d

, Bjørn Anders Pettersson Reif

^a

aDepartment of Technology Systems, University of Oslo, PO Box 70, Kjeller 2027, Norway

bNorwegian Defence Research Establishment, PO Box 25, Kjeller 2027, Norway

cNorwegian Defence Estates Agency, PO Box 405 Sentrum, Oslo 0103, Norway

dInstitude of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, Oslo 0315, Norway

a rt i c l e i nf o

Article history:

Received 11 January 2019 Revised 12 March 2019 Accepted 17 March 2019 Available online 22 March 2019 Keywords:

Shock-particle interaction Particle cloud

Particle-resolved simulation Pseudo-turbulent kinetic energy Volume averaging

a b s t ra c t

Weinvestigatetheshock-inducedflowthroughrandomparticlearraysusingparticle-resolvedLargeEddy SimulationsfordifferentincidentshockwaveMachnumbers,particlevolumefractionsandparticlesizes.

Weanalyzetrendsinmeanflowquantitiesandtheunresolvedtermsinthevolumeaveragedmomentum equation,aswevarythethreeparameters.Wefindthattheshockwaveattenuationandcertainmean flowtrendscanbepredictedbytheopacityoftheparticlecloud,whichisafunctionofparticlesizeand particlevolumefraction.WeshowthattheReynoldsstressfieldplaysanimportantroleinthemomen- tumbalanceattheparticlecloudedges,andthereforestronglyaffectsthereflectedshockwavestrength.

TheReynoldsstresswasfoundtobeinsensitivetoparticlesize,butstronglydependentonparticlevol- umefraction.Itisinbetteragreementwithresultsfromsimulationsofflowthroughparticlecloudsat fixedmeanslipReynoldsnumbersintheincompressibleregime,thanwithresultsfromothershockwave particlecloudstudies,whichhaveutilizedeitherinviscidortwo-dimensionalapproaches.Weproposean algebraicmodel forthe streamwiseReynoldsstress basedonthe observation thatthe separated flow regionsaretheprimarycontributionstotheReynoldsstress.

© 2019 The Authors. Published by Elsevier Ltd.

ThisisanopenaccessarticleundertheCCBY-NC-NDlicense.

(http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction

Interaction betweenshockwavesandparticlecloudsareofin- terestinanumberofdifferentnaturalphenomena, aswellasin- dustrialapplicationsandsafetymeasuressuchasshockwavemit- igationusingporousbarriers(Suzukietal.,2000;Chaudhurietal., 2013).Italsofindsapplicationsinheterogeneousexplosives(Zhang etal., 2006) andexplosive dissemination of powdersandliquids (Zhanget al., 2001;Milneet al.,2010; Rodriguez etal., 2017). In coalmines,enhanced orsecondaryexplosions duetocoaldust is amajor safetyconcern(Ugarte etal.,2017; Shimura andMatsuo, 2018).Shockwaveparticle interactionalsooccursinanumberof naturalphenomena, with volcanic eruptions (Bower and Woods, 1996) as the prime example. There are also astrophysical exam- plessuchasejectionofstellardustfromsupernovae (Silviaetal., 2012). Moregenerally, high-speedmultiphase flowhas important

∗ Corresponding author.

E-mail addresses: a.n.osnes@its.uio.no (A.N. Osnes), Magnus.Vartdal@ffi.no (M.

Vartdal), m.g.omang@astro.uio.no (M.G. Omang), b.a.p.reif@its.uio.no (B.A.P. Reif).

industrial applications, such as liquid andsolid fuel engines and fluidizedbeds.Gas-turbinesoperatinginregionswithsuspensions ofsand particlesin theair aresubjectto substantial degradation duetoparticledeposition onturbineblades (Hamedetal., 2006).

Waterinjectionsystemshavebeenusedtoreducesoundintensity atrocketlaunchpads(Ignatiusetal.,2008),anditmightbepossi- bletoutilizesimilarsystemstoreducejetnoise(Krothapallietal., 2003),whichisespeciallyimportantaroundair-craftsduringtake- off.

Shock wave interaction with particle clouds has been exten- sivelystudiedoverthelastfiftyyears.Thediluteparticlecloudand thegranularflowregimesarequitewellunderstood,buttheinter- mediateregimehasprovenchallengingtomodel(Theofanousand Chang, 2017). The intermediate volume fraction regime is where particlesneitherdisplaythesamestatisticalpropertiesasisolated particles, nor asparticles in thegranular regime. InCrowe etal.

(2011),flowswithparticlevolumefractionsabove 0.1%werecon- sideredto belong to thisregime, while Zhang etal. (2001) used a lower limit of1%. Thedifficulty in modeling theseflows stems from the complex interaction between the flow field and the https://doi.org/10.1016/j.ijmultiphaseflow.2019.03.010

0301-9322/© 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license. ( http://creativecommons.org/licenses/by-nc-nd/4.0/ )

particle distribution. The particles occupy a volume that is large enoughthattheircollectivenozzleeffectisoneofthedominating dynamical effects in the flow (Mehta et al., 2018b), but they are notcloseenoughthatparticlecollisionexclusivelydeterminesthe movement oftheparticles.Eachparticledeflectstheflow around it, causing local flow acceleration and deceleration that depends on the local particle configuration. Additionally, boundary layers develop over the particle surface, and the flow separates behind theparticle.When theshockwavepassesover aparticlethereis a reflection fromthe front of the particle, anda focusing ofthe shockwavebehindit(Tanno etal.,2003).The reflectedshockin- teractswiththeupstreamparticlesandtheirwakes,andalsowith reflected shocks fromother nearbyparticles. Thesecomplex flow dynamicsleadto alargevariation indragforcesthat dependson thelocalparticleconfiguration.

The intermediate particle volumefraction regime hasrecently become feasible to studyin much greater detail than wasprevi- ously possible.In experiments,theshorttime-scalesandthelim- ited possibility of recording data in the regions of interest have presented significant difficulties. Recent improvements to exper- imental techniques have enabled experimentalists to accurately characterize the wave system and particle distribution when a shock wave passes through a curtain ofparticles (Wagner et al., 2012; Ling et al., 2012). Even more recently, DeMauro et al.

(2017)usedhighspeed,time-resolved,particleimage velocimetry tomeasure velocityfieldsinfrontofandbehindtheparticlecur- tain.Thenewsetsofexperimentaldatahaveresultedinarenewed efforttostudytheseproblemsusingnumericalsimulations,inpar- ticularusingtheEulerian–Lagrangianframework(HouimandOran, 2016; Shallcross and Capecelatro, 2018; Theofanous and Chang, 2017), but also Eulerian–Eulerian models (McGrath et al., 2016;

Saureletal.,2017;Utkin,2017).

Some quantities are very difficult to measure experimentally, such asflow field distributionsinsidethe particlecloudandfluc- tuations at the particle scale. Instead, these can be obtained us- ing particle-resolved numericalsimulations. Such simulations are computationallyexpensive,sincealargenumberofparticlesmust be used to obtain meaningful statistics. A limiting factor is the very large scale separation between the dynamically important particlescalephysics andthegloballengthscaleoftheproblems.

However, anumberofstudies haverecentlybeenableto investi- gate shock-wave particlecloud interaction usingtwo-dimensional (Regeleetal.,2014;Hosseinzadeh-Niketal.,2018;Senetal.,2018) andthree-dimensional(Sridharanetal.,2015;Mehta etal., 2016;

2018a; 2018b; Theofanouset al., 2018) numerical simulations. In particular,particleresolvedsimulationscanbe usedtoinvestigate closures forunresolved termsthat appear inEulerian–Lagrangian orEulerian–Eulerian modelsduetoaveragingofproducts offluc- tuations.VolumeaveragingisoneformoffilterusedinLargeEddy Simulations,andisalsousedintheformulationofmostEulerian–

Lagrangian methods. Volume averaging does not commute with spatialortemporalderivatives.Therefore,the averagingoperation introduces newterms in the volume averaged equations, asdis- cussedine.g. SchwarzkopfandHorwitz(2015).Volumeaveraging canbeappliedtodatafromparticle-resolvedsimulationstoinves- tigate the unclosed terms, both the terms that appear in single- phasemodelsandthosethatarespecifictodispersedflowmodels.

Inthe intermediate particlevolumefraction regime,the inter- particle distanceandthe particlesize areof thesame order.The spatial extent ofthe flow field fluctuations is comparableto the inter-particle distance,andwewillrefer totheseasparticlescale fluctuations. It is common to divide the flow fluctuations into pseudo-turbulentandturbulentstructures,andrefertothekinetic energyintheseaspseudo-turbulentkineticenergy(PTKE)andtur- bulent kinetic energy (TKE), respectively. The flow perturbations inducedbytheparticlesareconsideredtobepseudo-turbulentef-

fects. Pseudo-turbulent flow structures might have very different timeandspatialscalesthanturbulentstructures.Forthisreasonit isausefultechniqueforanalysesandmodelingpurposestodistin- guishbetweenthetwo (Mehrabadi etal., 2015).In thesettingof shockwavespassing through particleclouds, PTKEis causedpri- marily by three effects. Firstly, shock wave reflection from indi- vidual particles causes very large differencesbetween theregion affectedbythe reflectedshockandthe surroundingregions,with correspondinglyhighPTKEvalues.Secondly,flowdeflectionaround particlescauses localflow accelerationsanddecelerations, result- inginbothstreamwiseandspanwisefluctuations.Inaddition,flow separationbehindparticlesalsocausesasignificantdeviationfrom the mean flow speed. This last effect in particular will be dis- cussed in this paper. As evident from these examples, pseudo- turbulent flow fluctuations are quite different fromclassical tur- bulence.However,pseudo-turbulentfluctuationsmightthemselves generateclassical turbulentfluctuations.Thisisexpectedtooccur asa resultofthestrongvelocity shearintheparticlewakes. The processesgeneratingTKEandPTKEareverydifferentphenomena, andshould thereforebe modeled differently.Both TKEand PTKE enterthe volumeaveraged momentum equationsthrougha term that is analogousto the classical Reynolds stress, andfor conve- niencewewillusethetermReynoldsstressforthistermthrough- outthispaper.

Thevelocityfluctuationsinshock-particleinteractionhavepre- viouslybeenexaminedintwo-dimensionalflowsusingbothinvis- cid(Regeleetal.,2014)andviscoussimulations(Hosseinzadeh-Nik et al., 2018). In those flows, the PTKE was found to be slightly higher in the inviscid simulations, but of the same order as the mean flow kinetic energy in both cases. Regele et al.(2014) ad- ditionally demonstrated the importance of capturing the PTKE in volume-averaged simulations in order to obtain correct pres- surefields. Incontrast, Mehta et al.(2018a) found very low val- ues of PTKE, demonstrating a significant difference between the two-dimensional and inviscid three-dimensional simulations. The Reynolds stress plays an important role in the dynamics around theparticlecloudedges,andinparticularitinfluencesthe(time- dependent)strengthofthereflectedshockwave.Sincetheincom- ing flowfield is alteredby thereflected shockwave, phenomena such as particle drag, pressure drop through the particle cloud andeven the transmitted shockstrength depend directly on the strengthofthereflectedshock,andthereforealsoontheReynolds stress.

Recent studies have recognized that an issue for Eulerian–

Lagrangianmethodsisthat theforcesimposedonthecontinuous phasebyaparticledisturbstheflowaroundtheparticle.Calcula- tionofdragbystandarddraglawsisincorrectifcontinuousphase variables within the disturbedflow region is used,because most drag-lawsarecalibratedagainstundisturbedflowquantities.Meth- odsforhandlingthisproblemexactlyinthezeroReynoldsnumber limit have been proposed andeven shown to yield good results infinite Reynolds numberflows (Horwitz andMani, 2016;2018;

Balachandar et al., 2019). Accounting for how the particle influ- encesthe localflow field wasalsodone in MooreandBalachan- dar(2018),who useda linearsuperposable waketo approximate continuousphase fluctuations within a particle cloud forincom- pressible flow. So far, no such model has been proposed for the compressibleflow inside a particlecloud. However, using knowl- edgeabouthowthe particlesdisturb the flowintheir vicinity to improvedragcomputationcanbedone eveninthiscomplexset- ting,aswillbeshown.

In this work we perform three-dimensional, time dependent, viscous simulations of a shock wave passing through a random particlearray. The particlesare assumed tobe inert andstation- ary.We varytheincident shockwaveMach numberbetween2.2 and3,theparticlevolumefractionbetween0.05and0.1,andthe

particlediameterbetween50μmand100μm.Preliminaryresults fromthesesimulationswerereportedinVartdalandOsnes(2018). Weutilizevolumeaveragingtodefinethemeanflowandthefluc- tuations from that mean.Key flow properties such as mean ve- locity,pressureanddensity,aswellastheReynoldsstress andits anisotropy,areexamined.

The purposeofthisworkistwofold.Firstly,weanalyzetrends in mean flow properties over different combinations of volume fractionsandMach numbers than those that havebeen reported previously.In addition, we vary theparticle diameter, which has notbeendonebefore.Thisanalysisimprovestheunderstandingof thebulk effectof particle cloud properties onshock wave parti- clecloudinteraction.Thesecondpurposeistoanalyzestatisticsof theReynoldsstressesandtheiranisotropy,aswellastheirimpor- tanceintheflowdynamics.Totheauthors’knowledge,thishasnot yetbeeninvestigatedforviscoussimulationsofthree-dimensional randomparticlearrays.Suchdataarecrucialtothedevelopmentof Reynoldsstressclosuremodelsforshockwaveparticlecloudinter- action.

This paperis organizedasfollows.InSection 2 we briefly in- troduce the basic flow patterns occurring when the shock wave passes over a group of particles. The governing equations and thevolumeaveraged equationsusedforanalysisare describedin Section3.Section 4describesthe computationalmethod andthe set-upoftheproblemsunderconsideration.Section5presentsthe simulationresults.Firstthegrid-qualityischeckedbyexamination ofparticle dragand resolution of viscous shear layers. Next, we discusstheshockwave attenuation asit passesthroughthe par- ticlelayer. We then examine the mean fields, i.e.density, veloc- ity,pressure andMach numberdistributions throughoutthe par- ticlelayer, anddiscusstrendsaswe changeparticlevolume frac- tions,particlediametersandincident shockwave Mach numbers.

Nextwediscussthevelocityfluctuationsandtheiranisotropy.We discussthe momentumbalance aroundtheupstream edgeofthe particlecloud, tohighlight thedynamicimportance ofthe veloc- ityfluctuations. The discussion ofthe simulation results is final- ized by an examination of the particle drag coefficients andthe average particle forces obtained in the different simulations. In Section6weutilizethedatafromtheresolvedsimulationstopro- posemodels that capturesomeof theobservedpropertiesofthe flow.Weprovideanalgebraicexpressionforcombinationsofpar- ticlediameterandparticlevolumefractionsthatresultinthesame shock wave attenuation. We also propose an algebraic Reynolds stress modelbased on theeffect of separatedflow behindparti- clesonvolumeaveragedequations,andcomparethismodeltothe streamwise velocity fluctuation intensity obtained inthe simula- tions.Finally,concludingremarksaregiveninSection7.

2. Shock-inducedflowaroundparticles

In thissection, we presenta brief overviewof the flow dur- ingandafter theshockwave passesovera group ofparticles. In Fig.1,numericalschlierenimages(Quirk,1997)andinstantaneous streamwise velocities are shown for a time-series in one of the simulations(caseVII)thatwillbedescribedinSection4.Thetime sequencecoverstheshockwavepatternandthesubsequentdevel- opmentofparticlewakes.Initially,a planarshockimpactsonthe firstparticle,andaregularshockreflectionisformed atthefront ofthe particle(first frame).As the shockpropagates, a Mach re- flectionisobtained, whichcan bediscernedinthe secondframe.

During this time, the pressure difference between the front and thebackoftheparticleisverylarge.Duetothepresenceofmul- tiple particles, the individual reflected shocks coalesce and form areflected shock(Boikoet al., 1997;Wagner etal., 2012), which over time becomes nearly planar and propagates upstream. The particlesalsocauseshockwavediffraction, asclearlyseenbythe

curvedfrontinthethirdframe.Behindeachparticle,theshockis focusedandahigh-pressureregioniscreated(fourthframe).

Viscousforcesbecomemoreimportantfortheflowaroundpar- ticles whenthe particleReynolds numberisreduced. Thisisrel- evant also for shock particle interaction, since smaller particles correspondto lowerparticleReynoldsnumbers.Henceforth,“par- ticle Reynolds number” and “Reynolds number” will be used in- terchangeably. Sun et al. (2005) showed that depending on the Reynoldsnumber,thehigh-pressureregioncreatedbyshockfocus- ingcancreatetemporarynegativedrag-coefficients.Thisphenom- enawasonlyobservedforparticleReynoldsnumbersoftheorder 10³.ForlowerReynoldsnumbers,viscousforcescounteractedthis effect,andthetotaldrag-coefficientremainedpositive.Astheyvar- iedtheparticleReynoldsnumberfrom4900to49,theimportance ofviscousforcesincreaseddrastically,andforthe lowerReynolds number,thelate-timeviscousforceswasalmosttwicethemagni- tudeofthepressureforces.

TheparticleReynoldsnumberisoftenusedtocharacterizethe flow,anditistypicallybasedonundisturbedflow,orincidentflow, quantities.Itislikelythatcharacterizationbasedonincident flow properties are lessreliable in the shock wave particle cloud set- ting,duetothegenerationofacollectivereflectedshockwave.The strength ofthisshockwave determines the propertiesofthe gas that enters theparticle cloud, and is highlydependent on prop- ertiesofthecloud. Inaddition,particles withinthecloud areex- posedtothepseudo-turbulentflowinducedbyupstreamparticles, whichisverylikelytoaffectstatisticalpropertiesofthewake.

Thedevelopmentofflowseparationandparticlewakes canbe seeninthemiddleframesinFig.1.Boundarylayers developover the particle surfaces, and the flow separates. The particle wakes arehighlydistortedduetothepresenceofotherparticles,evenfor thefirstparticlesinthecloud.ForisolatedparticleswithReynolds numbersintherange50−300,Nagataetal.(2016)foundthatfor increasing Reynoldsnumbers, the separation linemoves forward along the particle surface. They also found that higher Reynolds numbersresulted inlonger separatedflowregions. Withinapar- ticlecloud,thelength oftheseparatedflowregionissignificantly affected bythe presence ofother particles, asisthe caseforthe wakebehindtheleftmostparticleinFig.1.Thisphenomenacannot bedescribedsolelybasedonReynoldsnumberandMachnumbers.

Thebottom frameofFig.1showsa snapshotofthe flowover the whole particle cloud. The different effects discussed above are visible in this frame, occurring at different spatial locations throughouttheparticlecloud.

3. Governingequations

The gas-dynamic processes considered in this work are gov- ernedbytheconservationequationsofmass,momentumanden- ergy.Indifferentialformtheseare

∂

^t

ρ

+

∂

k

( ρ

^uk

)

=0, (1)

∂

^t

( ρ

^ui

)

⁺

∂

k

( ρ

^uiu_k

)

⁼⁻

∂

ip+

∂

j

σ

i j, (2)

∂

^t

( ρ

^E

)

⁺

∂

k

( ρ

^Euk+pu_k

)

⁼

∂

j

σ

i ju_i

−

∂

k

( λ∂

kT

)

^{, (3)}

where

ρ

^{isthe massdensity, uis the velocity, pis thepressure,}

σ

i j=

μ

(

∂

ju_i+

∂

iu_j−2

∂

ku_k

δ

i j/3) ^{is the viscous stress tensor,}

μ

^is

the dynamic viscosity, E=

ρ

^e+0.5

ρ

^uku_k is the total energy per unitvolume,eistheinternalenergyper unitmass,

λ

^isthether-

mal conductivity, and T is the temperature. We utilize the ideal gasequationofstate,with

γ

=1.4,andrelateinternalenergyand temperature by a constant specific heat capacity. We assume a powerlawdependenceofviscosityontemperature,withanexpo-

Fig. 1. Numerical schlieren images (top and bottom row) and streamwise velocities (middle row) in a cut plane, when a shock wave impacts on and passes through a cloud of particles. Flow direction is from left to right. Frames are taken at times (t −t 0)/τp= 0 . 13 , 0 . 69 , 1 . 25 , 1 . 79 , 2 . 86 , 4 . 81 , 6 . 77 , 8 . 73 , 10 . 63 , 12 . 47 , 26 . 82 . t 0denotes the time when the shock wave is at x = 0 (upstream particle cloud edge) and τpis defined in Eq. (7) . In the middle frames, the colormap is linear between −900 m/s (black) and 900 m/s (white). The red dashed square shows the location of the zoomed view in the upper rows. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

nentof0.76,andwerelatethethermaldiffusivitytotheviscosity byassumingaconstantPrandtlnumberof0.7.

Intheanalysisoftheresults,weconsiderthevolumeaveraged equations ofmotion. Theseare obtained by applyingthe volume averaging operator to Eqs.1, 2,3.We use thenotation ·for vol- umeaveraging,

·

^forphase-averaging, and˜·forFavre-averaging.

Thedeviationsfromthephase-averagedandFavre-averagedvalues are denoted by · and ·, respectively. Phaseaveraging and vol- ume averaging are relatedby

α

·

=·, where

α

^{denotes the gas}

phasevolumefraction.Weusethesymbol

α

^{pfortheparticlevol-}

umefraction.Theproblemunderconsiderationisstatisticallyho- mogeneousintheyandzdirectionsandthereforethevolumeav- eragedequationscanbe expressedinone dimension.Thevolume averagedequationsarethen

∂

^t

( α ρ )

⁺

∂

^x

( α ρ

^u˜1

)

^{=0, (4)}

∂

^t

( α ρ

^u˜1

)

⁺

∂

^x

( α ρ

^u˜1u˜1+

α

^p

)

=

∂

^x

( α σ

¹¹

)

⁻

∂

^x

α ρ

^R˜11

+1 V

S

pn₁dS−1 V

S

σ

1kn_kdS, (5)

∂

t

α ρ

^E˜

+

∂

x

α ρ

^E˜u˜1+

α

^p

^u˜1

=

∂

x

( α σ

11

^u˜1

)

⁻

∂

x

( α λ∂

xT

)

−

∂

^x

α ρ

^eu1

−

∂

^x

α ρ

^R˜11u˜1

+D^u+D^p+D^μ+D^ap+D^aμ. (6) Intheequationsabove,R˜₁₁=u₁u₁ isastressduetovelocityfluc- tuations,analogoustotheclassicalReynoldsstressandwereferto thistermasReynoldsstressthroughoutthispaper.Thecontinuous phaseboundaryisdenotedbyS,Vistheaveragingvolume,andthe integralsrepresenttheforcesactingontheparticlesurfaces.D^u=

−1/2

∂

^x(αρuiu

iu

1) is the turbulent diffusion, D^p=−

∂

^x(^αpu1) ^is the pressure diffusion, D^μ=

∂

x(α^u_jσ_j1) is the turbulent viscous diffusion,D^ap=

∂

^x(α^a1^p)isthepressure-diffusioneffectduetothe turbulent massfluxa=

ρ

^u

/

ρ

, andD^aμ=−

∂

^x(αa1σ11)is the

analogousviscous diffusion effect.An investigation ofthe energy balanceduringshockwaveparticlecloudinteractionisoutsidethe scopeofthiswork,butweincludetheequationforcompleteness.

TheReynoldsstressappearsinbothEqs.(5)and(6),andtheterms containingitrepresenttheforcesduetovelocityfluctuationsand theworkdonebythoseforces,respectively.Thosefluctuationscan be both shear turbulence and pseudo-turbulent fluctuations. The physicalprocessesrepresentedbytheReynoldsstresswillbe dis- cussedinthispaperinordertoguideclosuremodeling.

4. Computationalmethodandset-up 4.1. Computationalmethod

The simulations in this work are performed using the com- pressible flow solver “CharLES”, developed by Cascade Technolo- gies. The governing equations are solved with an entropy-stable schemeon a Voronoi-mesh (Breset al., 2018), and a third order Runge–Kutta method for time stepping. A discussion of entropy stableschemescanbefoundine.g.Tadmor(2003);Chandrashekar (2013).

4.2.Problemset-up

We perform numerical simulations of shock waves passing throughafixed cloudofparticles,withvaryingshockwaveMach number,particlesizeandparticlevolume fraction.Fig. 2showsa sketchofthe computationaldomainandthe particledistribution.

The particle cloud has length L, and spans the domain in the y andzdirections.WedenotetheparticlediameterbyDp.Thepar- ticleconfiguration inthe figure is the configurationused for the simulationswiththe largestparticlediameter. Thecomputational grid consists ofstructured grids around each particle, which ex- tend 0.2Dp out from the particle surface, and an approximately uniformVoronoi-gridintherestofthedomain.Withinthisstruc- turedregion,thecontrolvolumesizeincreasesgeometricallywith distanceto the particle surface. Fig.3 provides an impression of

Fig. 2. Sketch of the computational domain and the particle configuration used for the simulations with the largest particles. The particles are located at 0 ≤x ≤L , where L = 1 . 2 √ ³

4 mm, and the computational domain extends 2 L /3 upstream and L /3 downstream of the particle cloud. The axis directions are indicated at the ori- gin. The span-wise extent is set to a constant multiple of the particle diameter, so that y = z = 8 √ ³

4 D p, and therefore varies depending on the particle size.

Fig. 3. Illustration of the mesh around particle, where the faces of each control- volume in this cut-plane are shaded according to the direction of their normals.

There is a structured mesh around each particle extending 0.2 D pout from the par- ticle surface, and a Voronoi-grid in the rest of the domain. Note that the control volume sizes in this figure are adjusted for illustration purposes.

the mesh around each particle. The particle positions are drawn froma uniform random distribution,so that anyposition within 0≤x/L≤1hasequalprobabilityofcontainingaparticle.Weaccept aparticlepositionifitisnotcloserthan1.5Dp toanyotherparti- clecenter. Thisensures that the structured grids do not overlap, andthat there is a small distance between the structured grids where the Voronoi-grid can create a smooth transition between thetwostructuredregions.Inaddition,werequirethatstructured grids do not intersect the spanwise domainboundaries. Particles are drawn inthis wayuntil the particle volume fraction reaches the desired value. For the simulations considered in this study, theminimal numberofparticlesinanysimulationis586andthe maximal numberis 1173.The size of the control volumesin the Voronoipartof the meshmatches approximately theouter layer ofthestructuredgridaroundeachparticle,anditisslightlycoars- enedintheregions awayfromtheparticlecloud.The totalnum- ber of control volumesis roughly 6×10⁷ for all the simulations here.On300cores,eachsimulationtookroughly24hourstocom- plete. The initial condition consists oftwo homogeneous regions separated by a shock wave, where the pre-shock conditions are setto

ρ

⁰=1.2048kg/m³,u⁰=0 m/sand p⁰=1.01325×10⁵ Pa.

Table 1

The different simulations considered in this study and key parame- ters. Ma is the incident shock wave Mach number, Re p,ISis the par- ticle Reynolds number based on post incident shock values, and n is the number density. L sis the sight-length, as defined in Section 5.2 .

Case Ma L s/ L αp Re p,IS n [ mm ⁻³] D p[ μm]

I 2.2 0.196 0.1 6160 191.0 100

II 2.4 0.196 0.1 7091 191.0 100

III 2.6 0.196 0.1 7927 191.0 100

IV 2.8 0.196 0.1 8666 191.0 100

V 3.0 0.196 0.1 9309 191.0 100

VI 2.6 0.157 0.1 6292 382.0 79.4 VII 2.6 0.125 0.1 4994 763.9 63.0 VIII 2.6 0.163 0.075 4537 763.9 57.2 IX 2.6 0.226 0.05 3964 763.9 50

X 2.2 0.099 0.1 3080 1528 50

XI 2.6 0.099 0.1 3964 1528 50

XII 3.0 0.099 0.1 4654 1528 50

The post-shock conditions are determined from the shock wave strength.Post-shockquantitiesareusedfornormalization,andwill be denoted with a subscript IS. For the post-shock gas velocity we omit the numeric componentsubscript for notational conve- nience. The shock wave propagates in the x-direction. The up- streamboundaryissettothepost-shockconditionandthedown- stream boundaryis set toa zero-gradient outlet. We apply sym- metryconditionsattheyandzboundaries.

Table1providesanoverviewoftheparametercombinationswe simulate. The simulationswill be referred toasCase I,II, …,XII.

Wevary theincident shockwave Mach numberbetween2.2and 3,particlesizebetween50

μ

^mand100

μ

^{mandparticlevolume}

fractionbetween0.05and0.1.

The analysis is conducted using the volume averaging frame- work. We define averaging volumes spanningthe domain in the y and z directions, witha streamwise extent of L/60. These bins containboththegasphaseandtheparticles,andthustheparticle volumefractionwithinabinmightdeviatefromthebulkparticle volumefraction.Theflowquantitiesareaveragedoverthesebins, andwesubsequentlycomputeamovingaveragewithawindowof fivebinstoreducethesensitivityoftheresultstothelocalparticle configuration.

We utilizetwo time-scales to compare thesimulation results.

Theseare

τ

L=L

Ma

γ ρ

^p00

−1

,

τ

p=Dp

Ma

γ ρ

^p00

−1

, (7)

where

τ

L isthetimeittakesfortheincidentshockwavetotravel a distance equal to the particle cloud length,

τ

^{p is the time it}

takesfortheincident shockwave topass overa particleandMa istheMach number.Unlessotherwisespecified,thetime-scale is computed using the incident shock wave Mach number in each simulation, so that the time-scales are different for the differ- entsimulations. Welet t₀ denotethetime whenthe shockwave isatx=0.

5. Results

All simulations considered here feature the same basic flow pattern,whichhasbeenreportedinanumberofpreviousexperi- mentalandnumericalstudies.Themostimportantfeaturesarethe generationofthereflectedshock,thegenerationofparticlewakes, andthecontinuousweakeningoftheprimary shockasitimpacts onparticlesthroughoutthelayer.Itshouldbenotedthatthepar- ticlelayer considered hereisnot long enough tocompletely dis- sipatetheshockwave,andthereforeatransmittedshockemerges fromthe downstreamedge ofthe particlecloud. After thetrans-

−2 0 2 4 6 8 10 t/τp

0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

Cd

N = 492 N = 1002 N = 2252 N = 4842 N = 20252

Fig. 4. Drag coefficient for an isolated particle subjected to a Ma = 2 . 6 shock wave with different grid resolutions. Here, N denotes the number of points at the particle surface.

mittedshockhasmovedaway fromtheparticlecloud,aflow ex- pansionregionoccursaroundthedownstreamedgeoftheparticle cloud, where the flow transitionsfrom subsonic (Ma≈0.5−0.8) tosupersonic(Ma≈1.2−1.6).Theexpansionregionisterminated byashockwaveashortdistancedownstreamoftheparticlecloud.

5.1. Gridresolution

Previous studiesofshockinteraction withparticlearrayshave estimatedgridqualitiesbyexaminingthedrag-coefficientonasin- gle particle (Mehta et al., 2016; 2018b; Hosseinzadeh-Nik et al., 2018).Itisimportantthat theparticleforcesarewell reproduced inthesimulation,sincetheyarecentraltotheproblemunderin- vestigation.Followingthesameapproach,we conductsimulations ofasingleparticlewithdiameter63

μ

^m,subjectedtoaMa=2.6 shockwave withvarious numberoffaces at theparticle surface.

Thisparametercombinationischosen becauseitisinthemiddle of the rangeof Mach numbers andparticle sizeswe havesimu- lated. Fig. 4showsthe dragcoefficient, asdefinedinEq. (15), as a function of time for five different grid resolutions. The drag- coefficient with N=2252 deviates roughly 2% from the highest resolution, andisa feasible resolutionin termsof computational cost for the particle cloud simulations. Therefore, we apply this resolutiontothesimulationsconsideredinthispaper.

In addition to the drag-coefficient,we alsoexamine the grid- resolutionintermsoftheparameter

l⁺=

³

VCV/lviscous, (8)

whichisthenon-dimensionalgrid-lengthscalerelativetothelocal viscous lengthscaleoftheflow. Here

³

V_CV isthelengthscaleof thecontrolvolume,

lviscous=

μ

ρ (

^2Si jSi j

)

^{1/2, (9)}

is the viscous length scale (Wingstedt et al., 2017) and S_{i j}= 0.5(

∂

ju_i+

∂

iu_j)^{is the strain ratetensor. Theviscous length scale} can be interpretedasthe smallestlength scaleof thelocalshear flow. Thus, ifthe grid size iscomparable toor smaller than this length scale, it is reasonableto assume that the flow is well re- solved locally. The viscous length scale can be utilized indepen- dent ofthe state ofthe flow(turbulent orlaminar).It should be notedthatthevaluesobtainedforl_viscousdependsonthegrid-size,

0 20 40 60 80 100 120 140

l⁺ 0.000

0.005 0.010 0.015 0.020 0.025

Probabilitydensity

Fig. 5. Histogram of l ⁺ for case VII, with (Ma, αp, D p)= (2 . 6 , 0 . 1 , 63 μ^m), at (t −t 0)/τL= 1 . 36 .

andl⁺thereforeonlyservesasapost-simulationmeasureofgrid- quality,asopposed to avalue that can be used quantitatively to refinea mesh.Fig.5showsa histogramofthel⁺ valuesforcase VII,with(^Ma,

α

p, D_p)=(².6, 0.1, 63

μ

^m),at(^t−t₀)/

τ

L=1.36. Themiddle98%ofthedistributionislocatedbetweenl⁺=30and l⁺=114.Thehighestvaluesare locatedintheshearlayeraround eachparticleandtheirwakes.Forproblemswhereitiscriticalto resolvetheturbulentenergycascade,thegridsizeshouldbecom- parabletol_viscous,oriflargergridsizesareused,thesmallerscales should be appropriately modeled. In this problem, it is unlikely that details of the energy cascade are very important. Therefore therequirementonl⁺ canprobablybeslightlyrelaxedherewith- out affecting the results considerably. Some phenomena, such as wake-wakeinteraction,andshock-wakeinteraction,mightrequire finer resolutions than we use.However, it is not the purpose of this work to explore these phenomenain detail, andit is likely thattheyonlyhaveaminoreffectontheresultspresentedhere.

Theviscous length scales in thiscaseare distributedbetween 20nm and200 nm. Thismeans that thesmallestviscous length scalesareonlyaboutanorderlargerthantypicalmeanfreepaths ofairmolecules.TheKnudsennumberbasedontheviscouslength scaleisgivenby

Kn= l_free

lviscous = kBT

√2

π

^pDp,air

1

lviscous, (10)

wherel_free isthemeanfreepathofthemolecules,k_B istheBoltz- mannconstant andD_p,air=3.84×10⁻¹⁰mistheeffectivediame- ter ofan airmolecule.Throughoutmostof theparticlelayer, the Knudsennumbertakesvaluesaround0.1. Intheexpansion region atthedownstreamendoftheparticlecloud,theKnudsennumber increases,andaroundthe verylastparticles wefindvaluesup to two.Thisindicates thatthe wakesandshearlayers aroundparti- clesatthe downstreamendof theparticlecloud might beinflu- encedby non-continuumeffects, butwe do not expect those ef- fectstobeverylarge.

Weconcludethat thegrid-resolutionusedinthisstudyissuf- ficienttoobtainreliableparticleforces,andthat itrepresentsthe localflowgradientsinasatisfactorymanner.

5.2.Shockwaveattenuation

As theshockwave passesthrough theparticle cloud, it isat- tenuated by shock reflection from the particles. The amount of

0.0 0.2 0.4 0.6 0.8 1.0 x/L

−0.5

−0.4

−0.3

−0.2

−0.1 0.0

x/us−(ts−t0)[μs]

αp= 0.1, Dp= 100μm αp= 0.1, Dp= 79.4μm αp= 0.1, Dp= 63.0μm αp= 0.075, Dp= 57.2μm αp= 0.05, Dp= 50μm αp= 0.1, Dp= 50μm

Fig. 6. Difference between undisturbed shock arrival time, x / u s, where u sdenotes the incident shock wave speed, and the obtained shock arrival time, t s, as a function of position within the particle cloud for an incident Ma = 2 . 6 shock wave. The shock arrival time is defined as the time when the average pressure within the bin first exceeds 3 bar.

attenuationdependsonparticlesizeandparticlevolumefraction, as well as the regularity of the particle distribution. As will be shown,theshockwave attenuationandcertain meanflowtrends are well characterized by a single parameter depending only on geometricpropertiesoftheparticlecloud.

Fig. 6 showsthe shock arrival time as a function of distance withintheparticlecloudforthecaseswithMa=2.6.CasesVIand VIII,with (

α

^{p, Dp)} = (⁰.1, 79.4

μ

^m) ^and (⁰.075, 57.2

μ

^m) ^re- spectively, havevery similar shock speedsduring the passage of theshocksthroughtheparticleclouds.Cases IIIandIX,with(

α

p, Dp)=(⁰.1, 100

μ

^m)^and(⁰.05, 50

μ

^m)respectively,appearquite similarinthisplot,butitcanbeseenthatthedifferencebetween themincreaseswith distance.The shockwave reflectionimposes astrongtransientforceontheparticles,andthereforeshockwave attenuationservesasameasureoftheaverageinitialparticledrag, andviceversa.Thesimulations doindeedshowthatsummingup theforcesontheparticles duringthefirst

τ

p aftertheshockhits eachparticle,yieldsapproximatelythesameresultforcasesVIand VIII.

The regularity oftheparticle distributioncan affectthe shock attenuation through statistical differences in shock focusing and particleforces.Forthisreason,itisnecessarytoestimatetheeffect ofournon-randomparticle distribution.Wegeneratetheparticle distributions by drawing random positions satisfying two criteria that makes the particle distribution slightly more regular than a completelyrandomdistribution,asdiscussedinSection4.Toquan- tifythiseffect,we examine the area in they-z plane that isoc- cludedbytheparticlesasafunctionofdistance.Theoccludedarea atastreamwisecoordinatexistheprojectionofallparticleswith streamwise coordinates lessthan x,onto a plane, accounting for overlapbetweentheprojections.Werefertotheratioofthenon- occludedarea tothetotalarea asopacity,andcompute itforthe differentgeometries that we haveused inthe numerical simula- tions.Fig.7showstheopacityasafunctionofpositioninsidethe particlecloud, where each line is the meanof 8192 realizations.

The dashed lines are the corresponding opacities by only requir- ingthat the particles arecompletely inside theassigned domain.

The opacity for the less restrictive distribution is always slightly higherthan theopacity for ourparticle distribution.Our particle distributionsare possiblerealizationsofthelessrestrictive distri- bution, butas there is not too much overlap betweenthe stan-

0.0 0.2 0.4 0.6 0.8 1.0

x/L 0.00

0.25 0.50 0.75 1.00

Opacity

αp= 0.1, Dp= 100μm αp= 0.1, Dp= 79.4μm αp= 0.1, Dp= 63.0μm αp= 0.075, Dp= 57.2μm αp= 0.05, Dp= 50μm αp= 0.1, Dp= 50μm

Fig. 7. Opacity for the different geometry parameters used in the simulations, as a function of distance into the particle cloud. The solid lines are the results for the particle drawing method used in this study and the dashed lines are the results for distributions without the inter-particle distance restriction. The shaded areas indicate the standard deviations for each line. The dotted lines indicate L sfor the different parameter combinations.

darddeviationregions,itisclearthattheyareveryunlikely.Since ourdistribution occludesmore ofthe area overa givendistance, weexpect slightlystrongershockreflectionandaslightlyweaker transmitted shock.The effect of the restriction that the particles should not overlapthe spanwisedomain boundariesis examined in Section 6.1. The result is that it increases the opacity, and is moreimportantthantheinter-particledistance.Theresultswithin this work should be interpreted with these effects in mind. Ad- ditionally,itmustbe emphasizedthat onlya singlerealizationof the particle distribution is used forthe flow simulation for each parameter combination. Wealso note that the curveshaveslight bumpsnearx=0andx/L=1,duetotheconstraintsimposedon theparticledistribution.We expecta similareffectforthedistri- butionclosetothespanwisedomainboundaries.

The opacities for

α

^p=0.01, Dp=79.4

μ

^{m and}

α

^p=0.075, D_p=57.2

μ

^{mareverysimilar.Theseparameter}combinationsalso resulted in a very similar shock wave attenuation, which indi- cates that the opacity might be used to predict some proper- ties of shock-wave particle cloud interaction. Since the opacity curvesdo not seemto intersect,we use thelength at whichthe opacity equals 0.5asa unique numberthat represents theopac- ity. We refer to it as the sight-length, and use the symbol L_s. The sight-length for each configuration is given in Table 1, and markedinFig.7forthecaseswithMa=2.6.Weseethatthisclas- sification indicates that the parameter combinations (

α

^p, Dp)= (⁰.1, 79.4

μ

^m) ^and (

α

^p, Dp)=(⁰.075, 57.2

μ

^m) ^{should be} very similar, aswe do observe. However, theresults would indi- cate a larger difference between (

α

^p, Dp)=(⁰.1, 100

μ

^m) ^and (

α

p, D_p)=(⁰.05, 50

μ

^m) ^{than we observe. It should be noted} that there is considerable standard deviation in the sight length becauseofthenumberofparticleswe use,sotheapparent simi- laritybetweenthelattertwoparametercombinationscouldbeex- aggeratedbythespecificparticledistributions.Alargernumberof particlesoranensembleofsimulations couldbeusedtoexamine thisingreater detail,butthat isoutsidethe scopeofthecurrent work.

The results indicate that it is possible to characterize some propertiesofshockwaveparticlecloudinteractionusingthesight- length.Forthisreason,weprovideanalgebraicexpressionthatap- proximatesthisquantityinSection6.1.

0.0 0.5 1.0 x/L 0.4

0.6 0.8 1.0 1.2

u1/uIS

αp= 0.1, Dp= 100μm αp= 0.1, Dp= 79.4μm

αp= 0.1, Dp= 63.0μm αp= 0.075, Dp= 57.2μm

αp= 0.05, Dp= 50μm αp= 0.1, Dp= 50μm

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

Fig. 8. Mean flow velocity with Ma = 2 . 6 and different particle sizes and volume fractions at (t −t 0)/τL= 0 . 5 , 1 . 0 , 1 . 5 and 2.0 from left to right. The particle cloud is located between 0 ≤x / L ≤1.

0.0 0.5 1.0

x/L 0.4

0.6 0.8 1.0 1.2

u1/uIS

Ma= 2.2 Ma= 2.4

Ma= 2.6 Ma= 2.8

Ma= 3.0

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

Fig. 9. Mean flow velocity with D p= 100 μm , αp = 0 . 1 and different incident shock wave Mach numbers at (t −t 0)/τL= 2 .

5.3. Meanflow

Inthissection weexamine theflowfield duringandafterthe shock haspassed throughthe particle cloud. The flow quantities arephase-averaged,orFavre-averagedwhereappropriate,overvol- umesspanningthedomainintheyandzdirections.

Fig. 8 showsthe normalizedvelocity at (^t−t₀)/

τ

L=0.5, 1.0, 1.5and 2.0 forthe cases withMa=2.6. In the first two frames the shockwave islocated inside theparticlecloud. The reflected shock isvisible asa sharp dropin velocity slightly before x=0, andtherecoveryshockispresentaroundx/L=1.1inthelasttwo frames. For a given particle volume fraction, the reflected shock wave is strongerfor smaller particles. This is expectedbased on thebehavior oftheprimaryshockdiscussedabove, sincethereis ahigherattenuationoftheshockwaveinthesecases.Attheup- streamparticlecloudedge,themeanflowspeedincreasesrapidly over adistance equalto afew particle diametersandthen hasa gentler slope throughoutthe central regionof theparticle cloud.

At thedownstreamparticle cloudedge thereisa strongflow ex- pansion,andtheflowspeedroughlydoublesover0.9≤x/L≤1.05at latetimes.Asseeninthetworightmostframes,thestrengthofthe expansionincreaseswithtimeaftertheshockhasexitedtheparti- clelayer.Asnotedinthediscussionofgrid-sizeabove,theexpan- sionregionmightbesubjecttonon-continuumeffectsduetothe increasingmeanfreepathoftheairmoleculesovertheexpansion

region.Weobtainstrongerexpansionsforhighervolumefractions andsmaller particles. Thus it might be necessary to account for non-continuumeffects if the volume fraction isincreased or the particlediameterisdecreased.Theflowspeedvarieswithvolume fractionandparticlesizesinthesamemannerastheshockspeed discussed above. Again we find that the parameter combinations (

α

^p, Dp)=(⁰.1, 79.4

μ

^m)^and(

α

^p, Dp)=(⁰.075, 57.2

μ

^m)^are approximatelyequal,butthereareslightdifferencesbetweenthese casesintheregionupstreamoftheparticlecloud andwithin the expansionregion.

The variation of the mean velocity with incident shock wave Mach number is shown in Fig. 9. The normalized flow velocity withinthe particlelayerdecreases withincreasingMach number.

However,thenormalized velocitywithin afew particlediameters ofthe shock wave has a very weak dependenceon the incident shock Mach number. The expansion region accelerates the flow toabout1.2u_IS inall cases,andthereforethe relativestrength of the acceleration is larger for stronger incident shock waves.The reflectedshockhasa largerjumpinnormalized velocitywithin- creasingMachnumber.Wenote thatthestrengthofthereflected shock wave increases with time over the time-frame considered here(VartdalandOsnes,2018).

Itisofinteresttocomparetheresultsheretothoseobtainedin theinviscid simulationinMehta etal.(2018b).Fig. 10showsthe resultsfromcaseV,with(^Ma,

α

p, Dp)=(³.0, 0.1, 100

μ

^m),at

−9 −6 −3 0 3 6 9 12 15 x/Dp

0.0 0.2 0.4 0.6 0.8 1.0

u1/uIS

(t−t0)/τp= 4 (t−t0)/τp= 8 (t−t0)/τp= 12 Mehta et al. (2018)

Fig. 10. Comparison of normalized velocity for case V, with (Ma, αp, D p)= (3 . 0 , 0 . 1 , 100 μ^m), and the inviscid simulations of Mehta et al. (2018b) at three times. Black lines are the results from the simulations in this work, and orange lines are the inviscid simulation results. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

threetimesandthecorrespondingresultsfromtheinviscidsimu- lationsofMehtaetal.(2018b)withthesameincidentshockwave Machnumberandparticlevolumefraction.Onaqualitative level, thetwosimulationsagreequitewell.However,thereappearstobe anon-negligibledifferenceinthereflectedshockstrengthbetween the inviscid and viscous simulations, where the viscous simula- tionsfeatureastrongerreflectedshockwave.Weemphasizeagain thatourparticleconfigurationisslightlymoreregularthanthatin Mehtaetal.(2018b),whichweexpecttoaffectthereflectedshock strength.Thiswasalsoobservedwhencomparingthesimulations ofthe inviscidface-centered cubicarray inMehtaetal.(2016)to therandomarray.Wealsoexpectthattheviscouseffectsincrease thereflectedshockstrength,sincetherewillbebothstrongerpar- ticledragandalsoamuchstrongereffectofReynoldsstresses.The effectoftheReynoldsstresswillbediscussedbelow.

Themeanpressureanddensityprofilesare shownforthedif- ferentgeometriesinFig.11.Bothquantitiesdisplaymuchthesame behaviorasthevelocity field.There isarapidchange aroundthe upstreamparticlecloudedge,followedbyanapproximatelymono- tonic decrease throughout most of the particle cloud until quite close to the shockwave position. Forthe cases shown here, the pressure tends to the same level around the downstream parti- cle cloud edge, butthe pressure drop over the expansion region is significantly smaller for the lowest volume fraction case. The densityprofilesintersect insidethe particlelayer,andtheconfig- urationswithhigherLs havelower densities atlow xandhigher densitiesathigherx.Itcanbeseenthatthevelocity,pressureand densityprofiles vary predictably withthe sight-length. For lower Ls thereisa highermeanvelocity insidethe particlelayer, lower pressure and flatter density profiles. We find that case IX, with (^Ma,

α

^p, Dp)=(².6, 0.05, 50

μ

^m), deviates slightly from the trendobservedfortheothercases.Thismightindicatethatsome ofthetrendswe observemaybe slightlydifferentatlow volume fractions,oritmaybeaneffectofthespecificparticledistribution.

Fig.12 shows thelocal flow Mach numberforthe caseswith incidentshockwaveMachnumber2.6.Aswevarytheparticlevol- umefractionandparticlediameters,wefind,asexpected,thatthe Mach numberis lower forthecases wherewe observed a lower shockwave speedinsidetheparticlelayer.TheMachnumbersta- bilizes about 0.2L behindthe shockwave, andhas a slightposi- tivegradientovertheinteriorregionoftheparticlecloud.Thelo- calMachnumberdropstovaluesaround0.5–0.6forlate timesin these cases. However, it increases drastically over the expansion region,attaining valuesup to1.4inthe latestframe shownhere.

Thetransitiontosupersonicflowhappensaboutoneortwoparti- clediametersupstreamofthedownstreamcloudedge.

Fig.13 showsthe local flow Machnumber forthe caseswith

α

p=0.1, D_p=100

μ

^{m. When compared using the time-scale}

basedonincidentshockwave speed,wefindthathigherincident shock wave Mach numbers result in higher local Mach numbers within the particle cloud anddownstream. The expansion region is stronger forhigher Mach numbers,but itappears toconverge forthehighestvalues.However,whenwecomparethelocalMach numberatthesamephysicaltime, asshowninthelowerpanels,

0.0 0.5 1.0 1.5 2.0 2.5 3.0

p/pIS

0.0 0.5 1.0

x/L 0.5

1.0 1.5 2.0

ρ/ρIS

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

Fig. 11. Mean flow pressure (top) and density (bottom) with Ma = 2 . 6 and different particle sizes and volume fractions at (t −t 0)/τL= 0 . 5 , 1 . 0 , 1 . 5 and 2.0 from left to right. Line colors are as in Fig. 8 .

0.0 0.5 1.0 x/L

0.4 0.6 0.8 1.0 1.2 1.4 1.6

Ma

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

Fig. 12. Local Mach number for different geometries with incident shock wave Mach number 2.6 at (t −t 0)/τL = 0 . 5 , 1 . 0 , 1 . 5 and 2.0 from left to right. Line colors are as in Fig. 8 .

0.4 0.6 0.8 1.0 1.2 1.4 1.6

Ma

0.0 0.5 1.0

x/L 0.4

0.6 0.8 1.0 1.2 1.4 1.6

Ma

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

0.0 0.5 1.0

x/L

Fig. 13. Local Mach number with varying incident shock wave Mach numbers with αp= 0 . 1 and D p = 100 μm at (t −t 0)/τL = 0 . 5 , 1 . 0 , 1 . 5 and 2.0 from left to right (top).

The bottom panels show the results at the same physical time, corresponding to using the τLscaling with Ma = 2 . 6 , i.e. the same times as shown in Fig. 12 . Line colors are as in Fig. 9 .

we find that after the strongtransient following theshock wave passageateachlocation,all thesimulationsattainthesamelocal Mach number.Wealsoseethat thereflected shockwavehasthe same jumpinMach numberforthesesimulations. In theexpan- sionregion,theresultsdiffer,andtheincreaseinMachnumberis muchstrongerforhigherincidentshockwaveMachnumbers.

Regele et al. (2014) reported an average local Mach number about 0.4 for their inviscid two-dimensional simulations with a Ma=1.67 shock wave and a particle volume fraction of 0.15.

Our results indicate that the local Mach number does not de- pend much on the incident shock wave Mach number, but has a strong dependence on particle size. The lowest average Mach numberwithin thecloudinour caseswithMa=2.6happensfor

α

=0.1 and Dp=50

μ

^m, where the average value at late time is 0.55. We expect that, in addition to differencescaused by the two-dimensionality, theregularity oftheparticle configurationin

Regeleetal.(2014) strengthensthe reflected shockwave andre- sultsinalowerlocalMach numberthanfora randomconfigura- tion.

Insummary,wefindthatwithinthecentralpartoftheparticle cloud,variationofmeanflowpropertieswithparticlevolumefrac- tionandparticlediameteriswellrepresentedbythesight-length, whichcharacterizes thearea blockageper distance.Atthe down- streamedge,weobserveastrongflowexpansion,andwefindthat there is a region around the upstream cloud edge that behaves differently than the central region. In these regions, which have a streamwise extent of a few particle diameters,the behavior is not predicted by the sight-length. This is because the flow field fluctuationsaredynamically importantintheseregions,andthey dependdifferentlyonvolumefractionandparticlediametersthan the mean flow fields.This will be further discussed in the next section.