1886698

(1)

ContentslistsavailableatScienceDirect

International Journal of Multiphase Flow

journalhomepage:www.elsevier.com/locate/ijmulflow

Performance of drag force models for shock-accelerated ﬂow in dense particle suspensions

Andreas Nygård Osnes

^∗

, Magnus Vartdal

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

a rt i c l e i nf o

Article history:

Received 27 October 2020 Revised 17 December 2020 Accepted 6 January 2021 Available online 15 January 2021 Keywords:

Drag law Shock wave Dense gas-solid ﬂow Particle-resolved simulations Eulerian-Lagrangian simulations

a b s t r a c t

Modelsforprediction ofdragforceswithinaparticlecloudfollowingshock-acceleration areevaluated withtheaidofresultsfromparticle-resolvedsimulationsinordertoquantifyhowmuchthedisturbances introducedbytheproximityofnearbyparticlesaffectthedragforces.Thedragmodelsevaluatedhere consistofquasi-steadyforces, undisturbedflowforces, inviscidunsteady forces,and viscous unsteady forces.Twodenseparticlecurtaincorrectionschemestotheseforces,basedonvolumefractionandinput velocity,arealsoevaluated.Themodelsaretestedintwoways;firsttheyareevaluatedbasedonvolume- averagedflowfieldsfromparticle-resolvedsimulations;secondly,theyareappliedinEulerian-Lagrangian simulations,andtheresultsarecomparedtotheparticle-resolvedsimulations.

Theresultsshow thatbothcorrectionschemes significantlyimprovethe particleforcepredictions,but the averagetotalimpulseon theparticlesis still underpredictedby bothcorrectionschemes in both tests.Withthevolumeaveragedflowfieldsasinput,thevolumefractioncorrectiongivesthebestresults.

However,intheEulerian-Lagrangiansimulationsitisdemonstratedthatthevelocityfluctuationmodel, associatedwiththevelocitycorrectionscheme,iscrucialforobtainingaccuratepredictionsofthemean flowfields.

ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Simulations ofparticlemotionrequireaccuratemodels forthe interactionforcesbetweenparticlesandthesurroundingﬂuidﬂow.

While particle-resolved simulations give accurate predictions of drag andparticle movement, they are in mostcases too computationally demanding to be applicableto full-scale systemsof interest. Therefore, it is necessary to use less computationally ex- pensive methods, such as Eulerian-Eulerian (EE) simulations or Eulerian-Lagrangian (EL)simulations.The simpliﬁcationsmadeby thesemodelsintroducenewmodelingchallenges,sincethemeth- ods mustbe suppliedwithmodelsforthephysicalprocessesthat occuratscalessmallerthanthecomputationalgrid,suchasparti- cledrag. Forisolated particles,thereexistsa numberofdragcor- relations thatareveryaccurate,andcanmodeltheforcesreliably overalargerangeofMachnumbers,Reynoldsnumbersetc.How- ever,inthepresenceofotherparticles,itisquestionablehowwell these modelsare ableto capturethe particledrag, since they do

∗Corresponding author.

E-mail addresses: Andreas-Nygard.Osnes@ﬃ.no (A.N. Osnes), Magnus.Vartdal@ﬃ.no (M. Vartdal).

not account for the fluid mediated particle-particle forces.Addi- tionally, the classical dragmodels requirevalues forundisturbed fluidpropertiesattheparticlelocation,i.e.thefluidpropertiesthat would be observed there inabsence of the particle.These prop- ertiesare generallynot available inproblemswithdense particle suspensions.Therefore, it isnontrivial to compute particle forces in such problems, and while many studies have performed simulations of this kind, e.g. (Gai et al., 2020; Osnes et al., 2018;

Sugiyamaetal.,2019;Utkin,2017),itisunclearwhetherthephys- icalprocessesatparticlescaleare wellrepresented,andtherefore whetheror not theresults are reliable.This motivatesa detailed assessmentoftheperformance ofclassicaldraglawsinproblems featuringdenseparticlesuspensions.Thepurposeofthisworkisto assesshowwellclassicaldraglawscanrepresenttheparticledrag inthesettingofshock-acceleratedﬂowthroughalayerofstation- aryparticlesat10%volumefraction.

To characterize how well the drag laws can predict the drag forces,we usedatafromparticle-resolvedLargeEddySimulations ofa Mach 2.6shockwave propagating through a stationarypar- ticlelayer withan initial particle Reynolds numberof 2000. The datafromthesesimulationswere analyzedinOsnesetal.(2020). Inthatwork, itwasdemonstratedthatdirectapplicationofstan- darddrag-laws underpredictedthelate-timeparticleforcesbyup

https://doi.org/10.1016/j.ijmultiphaseﬂow.2021.103563

(2)

to 50%, and performed worse with increasing particle Reynolds number.Here,wetaketwodifferentapproachestoidentifytheap- plicabilityofthedraglaws.First,weusetheflow dataandparti- cleforcesfromtheparticle-resolvedsimulationsdirectly,andcom- parehowwelltheparticleforcesarepredictedbyclassicalmodels based onthe volume-averagedflow properties.Secondly,we per- formELsimulationsofthesameproblem,andagaincomparehow well theparticle forcesare represented.In thissecond approach, wealsocomparethevolume-averagedflowpropertiestothoseob- served intheparticle-resolvedsimulations.Sincethereisastrong two-way coupling between the particles and the fluid flow, the particleforcemodelhasastronginfluenceonthemeanfluidflow.

Particle-resolved simulations of flow through dense particle suspensionshaveproventobeavaluableapproachforunderstand- ingthephysicalprocessesoccurringintheinteriorofparticlesus- pensions.Suchsimulationsgive accesstofinelyresolvedtemporal and spatialdata inregions that are very challengingto probe in experimental studies. Particle-resolvedsimulations havetherefore becomepopularinrecentyears.Forshock-acceleratedflows,three- dimensional simulations were used in Mehta et al., 2016, 2018, 2019, 2020; Theofanous et al. (2018); Vartdal and Osnes (2018); Osnesetal.,2019,2020.Severaltwo-dimensionalparticle-resolved simulations have also been conducted, e.g. Regele et al. (2014); Hosseinzadeh-Niketal.(2018).

The particle-resolved studies have shownthat thelocal parti- cleconfigurationhas alarge influenceonthe particledragforces (Akikietal.,2016;Mehtaetal.,2018,2019).Thisisthecaseboth fortheinviscidshock-related forcesandforthelater quasi-steady dragforces.Theforcesimposedby theshockwavecan beampli- fied or weakened due to shock wave diffraction. Later, after development of particle wakes, there are strong velocity gradients oriented inall directions within the suspension. Therefore,parti- clescanbelocatedinregionsofbothveryhighandverylowflow speeds,withcorrespondinglystrongorweakdragforces.Anaddi- tionalcomplicatingfactoristhepossibilityoflocaltransonicflow regions, where larger voids in the particlesuspension can act as expanding nozzles that accelerate the flow fromsubsonic to su- personic in a limited region.Shocklet formation around particles insuchregionsisalsopossible.

Forces during the initial, primarily inviscid, shock-accelerated flow indense particle suspensions have beenanalyzed in Mehta etal.(2018,2019).Peakdragdecayswithdownstreamdistance,as wouldbeexpectedduetoshockwaveattenuation,butthemagni- tudewasalsofound todependdrasticallyonstatisticalproperties of the particleconfiguration; face-centered cubic orsimple cubic packingshadlargerpeakdragcoefficientsthanarandomdistribu- tion onaverage.Similar configurationeffects arealsolikely tobe thecaseforthequasi-steadydragcoefficients,althoughtotheau- thorsknowledgethishasnotyetbeeninvestigatedforhigh-speed flows.

High-speed flow through particle cloudshas alsobeen shown to be highly unsteady, and flow fluctuations can be significant in many configurations. These fluctuations consist of both pseudo-turbulent fluctuationsandclassical turbulent fluctuations.

While there isnot aclear-cut distinction betweenthese, pseudo- turbulent fluctuations are relatedto thedisturbance flow around a particle, which does not need to be turbulent, but still acts as a Reynolds-stress term under statistical averaging. The ki- netic energy contained in such fluctuations can be significant Regeleetal.(2014);VartdalandOsnes(2018);Osnesetal.(2019a); Mehta et al.(2020), and thisis important for modeling. For example,thepressureinELsimulationsareoftencomputedbysub- tracting the meankinetic energyfrom thetotal energy andsub- sequently employingthe equation of state, while the correctap- proachwouldbe tosubtractthetotalkinetic energy.Flow fluctuations also affectthe drag-forces, andrecent workshave charac-

terizedthedistributionofpeakforces,temporaldrag-ﬂuctuations, and the temporally averaged forces (Mehta et al., 2019; Osnes etal.,2020).

Thereareseveralpropertiesoftheparticlecloudthat haveef- fectsonthedragforces.Thesepropertiesinclude,butarenotlim- itedto,theparticlevolumefraction,thevolumefractiongradients, andthelocalparticle conﬁguration.Intotal,this complicatedde- pendence represents a formidable modeling challenge. However, ﬂuid-mediated particle-particle interaction models are currently actively researched and are able to capture some of these phe- nomena.NotablestudiesinthisdirectionareAkikietal.(2017b,a); Senetal.(2018);MooreandBalachandar(2019);Balachandaretal.

(2020).

Inthiswork,weaimtocharacterizemodelsthatarereadilyim- plementedincommonELcodesinusetoday,suchasthemodelby Parmaretal.(2010).Wealsocharacterizeaphysics-based correc- tiontotheinputtothedragmodels,whichwasrecentlyproposed byOsnesetal., 2019,2020.Fluid-mediatedparticle-particleinteraction models requiresignificantly more effortto implement,of- ten by relying on available particle resolved DNS data, andhave yet to be extended to the higher flow speeds of interest in this work. Thus they fall outsidethe scope of the currentwork. The dragmodelsconsideredherearealsoofinterestforEEsimulation strategies,where detailedparticleconfiguration datais notavail- able.A recentdiscussion of theproperties ofsuch simulations is foundin(Foxetal.,2020).

Improved simulationcapabilities forhigh-speed,dense, multi- phaseﬂowsareadvantageousforawiderangeofproblems.Afew examples are volcanic eruptions (Zwick and Balachandar, 2019), meteoroidbreakupMcMullanandCollins(2019),needle-freedrug delivery systems (Truong et al., 2006), solid and liquid fuel en- gines (Cai et al., 2003; Shimada et al., 2006; Bravo et al., 2015), explosion mitigation or explosive dispersal Zhang et al. (2001); Milneetal.(2010);Gottiparthietal.(2014),andnoiseattenuation onrocketlaunchpads(Ignatiusetal.,2008).

Thisarticleis structured asfollows.InSection 2, theparticle- resolved simulation data is used as input to the particle force models and the results are compared to the forces obtained in theparticle-resolved simulations.Section 2.1brieﬂy describesthe particle-resolved simulations, Section 2.2 introduces the different forcemodels.InSection2.3,particle-resolvedsimulationsandthe force models are applied to two single-sphere problems, while Section 2.4compares the forcesof the shock-wave particle-cloud simulationswiththeforcepredictions.Next,theforcemodels are evaluatedaspartofELsimulationsinSection3.TheELsimulation method is presented in Section 3.1,while the simulation results andcomparisontotheparticle-resolvedsimulationsarecontained inSection 3.2.Finally,Section 4containsconcludingremarksand adiscussionofpossibleimprovementsoftheforcemodels.

2. Forceestimationfromparticle-resolvedﬂuiddata 2.1. Particle-resolvedsimulations

The particle-resolved simulations that are used as reference data in this work were described and analyzed in Osnes et al., 2020.Forthereadersconvenience, werepeatthemostimportant detailshere.Theproblemconsideredisthatofthepropagationof aninitiallyMach2.6shockwavethrougharandomarrayofparti- clesat10%volumefraction.Thesimulationsarethree-dimensional andincludeviscousterms,andtheparticlesareassumedtobesta- tionaryandinert.TheparticleReynoldsnumberis

Rep,IS=

ρ

ÎSûÎS^D^p

μ

^IS ⁼²⁰⁰⁰^, ⁽¹⁾

(3)

Fig. 1. Sketch of the problem setup and the computational domain. A particle layer with particle volume fraction 0.1 is located between 0 ≤x ≤L . A shock wave with Mach number 2.6 is initially located at −0 . 1 L and propagating towards the particle layer. The particle Reynolds number based on post-incident shock values is 20 0 0.

where subscript IS denotes a ﬂow variable behind the incident shock wave,

ρ

îs ^the ^fluid ^density, û îs ^the fluid velocityin the downstream direction, Dp is the particle diameter, and

μ

^is ^the

ﬂuidviscosity.Fig.1showsasketchoftheproblemsetup.Thepar- ticlelayerislocatedbetween0<x<L,wherexisthestreamwise coordinate and L is the particle layer length, which is 12√³

16Dp. Thespanwiseextentsofthedomainare8√³

4D_p.Thespanwisedo- mainboundariesareperiodic,whiletheupstreamboundaryisset to the post-shock state andthe downstream boundaryis a zero- gradient outlet. Ten simulations with different random distribution realizations were performed, and the results were volume- averaged as well asaveraged over this simulation ensemble.The initialpositionoftheshockwavewassetto−0.1L,andeachsim- ulationwasrununtilt=3.75

τ

L,where

τ

^L⁼ _Ma^L

ISc0, (2)

is thetimeit takesforthe initialshockwave totravela distance equal to the particle layer length. Here, c₀ isthe ambient speed of sound. The ﬂuid equation of state is set to an ideal gas with adiabaticindex

γ

=1.4.

The computational mesh consistsof an unstructured Voronoi- based grid with a control volume length scale of approximately 0.036Dp for regions closer than Dp/2 to any particle. In the rest ofthedomain,thecontrolvolumelength scaleisdoubled,except from part of the upstream anddownstream regions, where it is doubled again. The total number of control volumes is approxi- mately10⁸.Thisresolutionwasabletoaccuratelycapturethedrag onan isolatedparticleandtheviscouslengthscales werereason- ably resolved at Re_p_,_IS≈5000 (Osnes et al., 2019a). The simulations werealsoclosetoconvergedintermsofthevelocity fluctuation levels (Osnes etal., 2019b). Sincethe currentstudyconsid- ers Re_p_,_IS=2000, the grid-requirements are lessstrict, andthus theflowfield isbetterresolvedinthecurrentstudythanthatfor whichthegrid-resolutionwastested.Additionally,weconducttwo single-spheresimulationshereinordertoverifytheabilityofthe particle-resolved simulations to capturethe particleforces.These simulations arepresentedinSection2.3.Forfurtherdetails about thegoverningequations,thecomputationalmethod,andotherde- tails aboutthe simulations,the readerisreferred to Osnesetal., 2020.

The particle-resolved simulationdata containstheforce histo- riesfor9310particles.Theflowfielddatawasrecordedbyvolume- averagingoverbinswithastreamwiseextentofL/60≈0.5Dp and spanning the y andz directions.The averaged flow field is thus onlyafunctionofthestreamwisespatialcoordinatexandtime.All terms in the volume-averaged flow equations (see Osnes (2019)) were stored. The volume-averaged results thus contain the data that wouldbeavailable inan ELorEE simulationwitha stream- wisegrid-spacingequaltothebin-width.Inaddition,theparticle- resolved results include terms such as the correlations between

Fig. 2. Example force histories at different locations. These are taken from the particles in the ensemble that have streamwise coordinates closest to the six locations x/L = 0 , 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8 , 1 .

velocityfluctuations,thecorrelationbetweenpressure-fluctuations andvelocityfluctuation,etc.,that cannot bedirectlycomputedin an ELsimulation. These additional terms give information about the disturbance fields induced by the particles. In Osnes et al.

(2019,2020),asimplemodelfortheaverage velocitydisturbance field(orvelocity componentofthepseudo-turbulentfluctuations) was proposed. This model accounts for the fluctuations that are duetoparticlewakes, andwascalibratedusingthevelocity fluc- tuationcorrelations.Thismodelwillbeusedbelowtoapproximate theundisturbedflowfield,whichisneededinthedragforcemod- els.

2.2. Particleforcemodels

Examplesoftheforcehistoriesthat aretobeapproximatedby thedragmodelisshowninFig.2.Theﬁgureshowsthedragcoef- ﬁcientinthestreamwisedirectionforthesixparticlesintheen- semble that are theclosest to thesixstreamwise locationsx/L= 0, 0.2, 0.4, 0.6, 0.8, and1,asa function oftime. Here, time is normalized by the time it takes for the incident shock wave to traveloneparticlediameter,labeled

τ

p andgivenby

τ

^p⁼_Ma^D^p

ISc0, (3)

whilethedragcoeﬃcientisgivenby C_D_,_IS=

Sp

pn1−

σ

1jn_j

dS 0.5

ρ

^IS^u²ISπ

4D²_p , (4)

(4)

whereSp istheparticlesurface, pisthepressure,nisthenormal vectorattheparticlesurface,repeatedsubscriptsimplysummation over components 1–3, where component 1 is in the streamwise direction, anddSisasurfaceelement.Foreachparticle,theforce historystartswithaverysharppeakwhentheshockwaveimpacts ontheparticle.Subsequently,theforcedecaysrapidlyastheshock wave movesfurtherdownstream.After about10−20

τ

p,theforce attainsamorestablevalue,buttherearestillsigniﬁcantvariations abouttheslowlyvaryingmeanforce.

In Parmar et al. (2012), the following force decomposition is usedtosplittheforcesintodifferentphysicaleffects,

Fp=Fqs+Fun+Fiu+Fvu, (5) where F_p is the total force acting on the particle, while the four componentsare thequasi-steadydrag, theundisturbed ﬂuid (stressdivergence) force,the inviscidunsteady force,andthevis- cous unsteady force. Annamalai and Balachandar (2017) derived theexpressions forthesebasedonageneralizedFaxén’stheorem.

Theseare

Fqs=3Dp

πν ( ρ

^u

)

^un^S

(

^Re^p^,^Ma^p^,

α

^p

)

^, ⁽⁶⁾

where

ν

îs^the^kinematic^viscosityôf^the^fluid,ôverline^denoteâv-

eragingovertheparticlesurface(superscriptS)orvolume(super- scriptV),superscript“un” denotesundisturbedﬂuidproperties, isadragcorrelationfactor,typicallygivenbyanempiricalcorrela- tion,andMap=(^u^V−up)/c^V istheparticleMach number,where u_p is theparticle velocity andc isthe local speedof sound, and

α

^p îs ^the ^particle^volume ^fraction.Quasi-steady drag correlations areoftengivenintermsofthedragcoefficient,whichisafunction offar-fieldflowproperties.Inthecaseofauniformflowfield,the relationbetweenând^the^drag^coefficientîs

=Rep

24C_D_,_qs. (7)

Theundisturbedﬂuidforceisgivenby Fun=Vp

ρ

^un^Du_Dt^un

V

, (8)

whereVp istheparticlevolumeandD/Dt isthematerialderivative basedontheundisturbedﬂuidvelocity.Theinviscidunsteadyforce is

F_iu=Vp

t

ξ=−∞K_iu

(

^t⁻

ξ

^,^Ma^p

) (

^t⁻

ξ )

D

Dt

( ρ

^u^r

)

^un^S

t=ξ

d

ξ

^, ⁽⁹⁾

where K_iu is the inviscid unsteady kernel, while the viscous un- steadyforcesare

Fvu =18R²_p√

πν

^t

ξ=−∞K^V_vu

(

^t⁻

ξ

,Rep,Map

)

_dt^d

( ρ

^ur

)

^un^S

t=ξ

d

ξ

+6R²_p√

πν

^t

ξ=−∞K_vu^S

(

^t−

ξ

,Rep,Map

)

_dt^d

( ρ

^u

)

^un^S

t=ξ

d

ξ

,

(10) where K_vu^S is the viscous unsteady kernel for the surface contribution of theviscous unsteady force, K_vu^V isthe kernelrelatedto thevolumecontribution,uristheradialcomponentthevelocityas observed bytheparticle,and _dt^d isthematerialderivative following the particle.In these expressions, the ﬂow properties in Re_p andMap shouldbetakenastheaveragesoftheundisturbedﬂuid propertiesovertheparticlevolume.

We use the following models for the different force components. For the quasi-steady drag, we use the model by Parmaretal.(2010),where

CD,qs(^Rep,Map)

= C_D_,_std(^Rep)+[CD,Mcr(^Rep)−C_D_,_std(^Rep)^]^MaMcr^p if Map≤Mcr

C_D_,_sub(^Rep,Ma_p) ^if ^Mcr≤Ma_p≤1 CD,sup(^Rep,Map) ^if ¹<Map≤1.75

. (11)

Thereaderis referredto Parmaretal.(2009)fortheexpressions forthevariousparameters inEq.(11).Theundisturbedﬂuidforce is predominantly the pressure gradient force, and will thus be computedassuch.However,itisimportanttonoticethatitisthe undisturbedﬂowpropertiesthatenterinthisforceaswell.

Theinviscidunsteadyforce ismodeled withtheMach-number dependentkernels obtainedby Parmaretal.(2008,2009) intab- ulated form. It should be noted that the tabulated kernels are based on constant acceleration of a sphere in a fluid at a given backgroundMachnumber,wherethefluidaccelerationmagnitude is significantly less than that of the incident high Mach number shockwave considered here. Therefore,itis not obviousthat the kernelcaptures the relevantflow physics inthe presentconfigu- ration.Indeed,theresultsofParmaretal.(2009)indicatethatthe accuracyoftheforcemodelishighlydependentontheshockwave Machnumber(seeFig.6inParmaretal.(2009)).

Fortheviscousunsteadykernels,weusethemodelofMeiand Adrian(1992),wherethesurfaceandvolumecontributionsarenot givenseparately,butratherasinglekernelisused.Theviscousun- steadyforceisthen

Fvu=3

πμ

^D^p^t

ξ=−∞Kvu

(

^t^,

ξ )

_dt^d

( ρ

^u

)

^un^S^d

ξ

^, ⁽¹²⁾

where

Kvu(^t,ξ)=

4π(^t−ξ)ν D²_p

1/4

+

π|^u(ξ)−up(ξ)|³ Dpν^fH³(^Rep) (^t−ξ)²

1/2−2

,

(13) where f_H=0.75+0.105Re_p.Likefortheinviscidunsteadykernel, itshouldbenotedthattheviscousunsteadykernelisobtainedin adifferentﬂow regimethan appropriate forthe currentproblem, withsmalloscillationsoftheinﬂow.

Theaboveforcemodelsarederivedforisolatedparticles,where theparticlevolumefractionisnegligibleandtheundisturbedﬂow velocity can be easily estimated. Neither of thesestatements are true for the present case. To account for these differences, two model corrections are considered. For some of the results presented below, we will use the model proposed in Osnes et al.

(2019,2020)to approximate an undisturbedﬂowvelocity ateach particlelocation. The model introduces a correction factorto the volume-averagedvelocity,whichisafunctionofparticleReynolds numberandparticlevolumefraction.Themodelapproximatesthe undisturbedstreamwiseﬂowvelocityby

u^un=u˜

α

⁻

α

^sep

( α

^p^,^Re^p

)

^, ⁽¹⁴⁾

^sep ûsing â^history întegral ôver

the relative ﬂow velocity for all particles in the control volume.

However, such a model hasnot yet beendeveloped, and is out- side thescope of thecurrent work.Along withthe correction to

R=1.2635kg/m³, u_R=0, p_R=102203Pa, (19)

where

ρ

L, u_L, p_L denotethedensity,velocityandpressureonthe left side ofthediscontinuity, while

ρ

R, u_R, p_R denotethe corre- spondingvaluesattherightside ofthediscontinuity.Thediscon- tinuityislocatedatx=1.25D_p,andtheparticleatx=0.

Fortheshock-particleinteraction,theinitialconditionconsists ofthestates

ρ

^L⁼¹^.⁶⁵⁸²^kg^/^m³^, ^u^L⁼¹¹⁴^.⁴⁷^, ^p^L⁼¹⁵⁹⁰⁶⁰ ^Pa^, ⁽²⁰⁾

ρ

R=1.2048kg/m³, uR=0, pR=101325Pa, (21) and the particle diameter and gas viscosity are set by requiring Rep=4900.Forbothcases,resultsforgridswithsimilarresolution to that used in the full particle cloud simulations are presented.

Grid convergence studies for similar conﬁgurations are found in Osnes(2019).

The simulation results are shown in Fig. 3 along with the (digitized) results of Annamalai and Balachandar (2017) and Sunetal. (2005),andthe force modelpredictions. Theforce obtained in the particle-resolved simulation of the expansion fan agreesperfectlywiththeDNSresultsofAnnamalaiandBalachan- dar (2017). This is also the case for the force model prediction, which, in thiscase, is based on the exact solution ofthe corre- spondingshock-tubeproblemwithouttheparticle.

In the shock-particle case, the agreement of the cur- rent particle-resolved simulation and the simulations of Sun et al. (2005) is excellent up to 2tc_IS/Dp=6. At later times, the currentsimulation predicts a slightlyhigher force, whichwe attribute to the development of the particle wake, which can be expected to differ for axisymmetric and three-dimensional simulations. Compared to the experiments, slightly higher drag is obtained in the simulation. This is to be expected since the experiments were conducted at Rep≈3×10⁵, and thus the vis- cousforcecontributionisnegligibleintheexperiments.Theforce modelpredictionisalsoingoodagreement withthe simulations, althoughnottothesameextentasintheinviscidexpansioncase.

Nevertheless,consideringthattheforce comparisonspresentedin Parmar et al. (2009) for shock-particle interactions also showed differences between the model prediction and the simulations of Sun et al. (2005), and that the viscous unsteady kernel is developedforacompletelydifferentﬂowregime,weconsiderthe forcemodelpredictionshowninFig.3tobequitegood.

Inconclusion,the particle-resolvedsimulations areable toac- curatelycapturetheforcesonisolatedparticles.Theforce models alsocapturetheseforcesverywell.Theresultsforforcesindense particlesuspensionsarepresentednext.

2.4. Comparisonofforcemodelsandparticleresolvedsimulation results

Individual particle histories are interesting to examine since thesegiveanimpressionaboutthevariedbehaviourthatthedrag modelsshouldideallypredict.Fig.4showsthreeparticleforcehis- toriesalong withthedragforce predictions basedonthe models described above. Each particle history is shifted by a time t₀_,_p, which is the time when the force on that particle first exceeds 1‰ofthemaximalparticleforcemagnitude.Theparticlessharea commongeneralbehaviour,withasharpspikeintheforce,corre- spondingtotheshock-particleinteraction,followedbyarapidde- cayandthenamoregradualslowdecayovertime.Therapiddecay partisverydifferentforthesethreeparticles,whereespeciallythe particle shown in the bottom panel has a behaviour that is sig- nificantlyaffected by fluid-mediated particle-particleinteractions.

Oscillationswithatime-scale ofabout2

τ

^p−4

τ

^p ^are^clearly^seen,

andtheir amplitudes take values up to 50% of the average drag, seee.g.theforcearound(^t−t₀_,_p)=14

τ

^p^for^the^particle^shownⁱⁿ

theupperpanels.

(6)

Fig. 3. Particle forces for (a): an isolated particle in a weak inviscid expansion fan, (b): an isolated particle exposed to a Ma = 1 . 22 shock wave at Re p= 4900 . Here, c IS

denotes the speed of sound based on the post-shock state.

Fig. 4. Three particle histories as well as the drag force predictions. Solid black line: particle-resolved force, solid grey line: F qs+ Fun+ F iu+ F vu, dashed black line: F qs, dashed grey line: Fun, dotted black line: Fiu, dotted grey line: F vu. (a): isolated particle model, (b): velocity corrected model, (c): volume fraction corrected model.

The initialspikeintheforce imposedbythe shockwave, visi- bleinthehistories,ispredictedbythedragmodels.Itsmagnitude is overestimatedforall three particles withall three approaches, exceptfortheparticleforceshowninthetoppanelwiththeiso- latedparticlemodel.Thevolumefractioncorrectedmodelpredicts uptotwicethemagnitudeobservedintheparticleresolvedsimu- lations,whilethetwootherapproacheshavemoremodestoveres- timations.Theprimaryfactorsresponsible forthehigherpeaksin the volume-fractioncorrected model is the strong ampliﬁcations ofFqs andF_iu inthevolumefractioncorrected model,since these peak ataroundthesametime.Thisisinlinewiththe comments

madeabove,andconsistentwiththebehaviourobservedinFig.6b of(Parmaretal., 2009)ofscaling issuesofthekernelwithMach number.

Some particle force histories are better predictedthan others.

Forexample,theparticleforcepredictionintheupper-mostpanel fortheisolatedparticlemodelﬁtsbetterwiththeparticle-resolved force than the one in the middlepanel. This is expecteddue to the variation related to the particle conﬁguration. This variation is not modeled with the draglaws evaluated here, and compar- ison ofindividual particle histories is thereforenot a well-suited approach forevaluating the applicability of the draglaws to the

(7)

Fig. 5. Average forces. Black lines show the results for the particle-resolved simulations and the grey shaded regions show one standard deviation to each side. Colored lines show the average force model predictions, and the colored shaded regions show one standard deviation to each side. (a): Isolated particle model, (b): velocity corrected model, (c): volume fraction corrected model.

Fig. 6. Average force components. Black line: Fqs, blue line: F un, orange line: F iu, green line: Fvu. The shaded areas indicate one standard deviation. (a): Isolated particle model, (b): velocity corrected model, (c): volume fraction corrected model. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

currentproblem.Forthisreason,wealsoconsidertheaveragepre- dictionsfortheentirecloud.

Overall, none ofthe three approachesare able to capturethe average drag forces very well. The most signiﬁcant deﬁciencies are foundintheaverageforce magnitudes.Fig.5showstheforce history averaged over all particles, where the time is shifted by t₀_,_pforeach particle,forthe particle resolvedsimulationsandthe threemodelingapproaches.Qualitatively,theinitialspikeandsub- sequentdecayarecapturedbyallthreeapproaches.Thepeakforce is too high on average for the velocity corrected and the vol-

ume fraction corrected models, while it is too low for the isolatedparticle model. The peak force occurs slightlytoo late, but this is most likely related to the grid size used in this assess- ment.Atlatertimes,thepredictedmeanvaluesarelowerthanthe particle-resolved meanvaluesforall threeapproaches.The mean dragisapproximatelyatonestandarddeviationbelowtheparticle- resolvedmeanforcewiththeisolatedparticlemodel.Withtheve- locitycorrectedmodel,thisdifferenceis morethanhalved,while thevolumefractioncorrectedmodelisbetteragain.

(8)

The force variationspredictedby the dragmodels aresmaller thantheone fromtheparticle-resolvedsimulations.Thevariation islargestforthevolume-fractioncorrectedmodelandsmallestfor theisolated particlemodel.Underpredictionoftheforcevariation isexpected,becausethedragmodelwillpredictapproximatelythe sameforceonallparticles thathavesimilarstreamwisepositions.

The particle-resolved simulations alsohave a force variation that isrelatedtothelocalparticleconﬁguration.Consideringtheentire particle cloud, slightlymore than half thevariation of the forces atmosttimescan be explainedby thedifferencesindragatdif- ferentpositions,andslightlylessthanhalfofthevariationbythe particle-conﬁgurationrelateddragdifferences.

Thevarious forcecomponentsasafunction oftime areshown inFig.6.Theforce hasnon-negligiblecontributionsfromallcom- ponents inthedecomposition,Eq.(5),buttheunsteadyforcesare primarilyimportantinashortintervalaftertheshockimpact.Over thetimet−t₀_,_p∈[0,

τ

^p^],theinviscidforces(FunandF_iu)dominate.

Themagnitudeoftheinviscidunsteadyforceissignificantlylarger thantheothercomponentsduringthisinterval.Itshouldbenoted that theforcesthatinvolvegradients (F_un andF_iu)willbe slightly smoothedduetothecoarsespatialsamplingoftheflowfieldrel- ativetothethicknessoftheshockwave.Theviscouscontributions begintobeimportantataroundt−t₀_,_p=

τ

p.Theviscousunsteady force has the smallest magnitude overall. The quasi-steady force dominatesthedragaftert−t₀_,_p=2

τ

^p^.

There are clear differences in the relative importance of the force componentsforthethreeapproaches, andthescalingofthe inviscidunsteadyforceappearstobeinappropriate.Inthevolume- fractioncorrectedmodel,Fqs,F_iu andFvuaresimplyscaleduprela- tive totheisolatedparticlemodel,givingthemlargermagnitudes relative to Fun. First of all, considering that the peak drag force is higher than the particle-resolved results, this scaling appears to be inappropriate for the period of time following shortly after the shockwave passesover each particle.This isnot entirely surprising, since the conﬁguration information is unavailable on thetimescaleoftheshockinteraction,andtheﬂowconditionsare vastly differentthan forthose wherethe correctionsare derived.

Secondly, at the later stages of the process, the dragpredictions aretoolowandthusascalingofpositiveforcesappearstobeap- propriate.However, atthisstage F_iuisnegative, andthereforethe scaling ofthisforceactsina waythat brings thedragprediction furtherawayfromthelevelobservedintheparticle-resolvedsimulations. Thevelocitycorrectionprimarily affectsthequasi-steady dragforce,butitalsohasasmalleffectontheshapeoftheinvis- cidunsteadyforcehistory.Likethevolumefractioncorrection,the velocity correction modelwasobtainedforthe quasisteadycase, anddirectapplicationtotheunsteadyhighlytransientcomponents doesnotappeartobeappropriate.

The force models undershoot the averagedrag ataround (^t− t₀_,_p)≈2

τ

^p−6

τ

^p, whichis duetoF_iu.The inviscid unsteadyforce model produces a negative force contribution after (^t−t₀_,_p)≈3, resultingin underpredictionofthetotal force.Thisnegativeforce ispredictedbythe modelbecauseit isbasedontheﬂowaround an isolatedparticle,wherethenegativeforce isaresultofshock- wave diffractionandgeneration ofa high-pressureregion behind theparticle(Sunetal.,2005).However,inarandomparticlearray, theparticlepositionsareoftensuchthatmanyparticleshaveother particles in their immediate proximity, which affects the shock- particleinteractionthroughﬂuid-mediatedparticle-particleforces.

Thiscanpotentiallycanceloutthenegativeinviscidforces.Onav- erage,itappears thatthenegativeinviscid forcesaresigniﬁcantly dampened by the particle distribution, since no local minimum is found in the average particle-resolved simulation results.This effect is likely to be a function of the average inter-particle distance and/or volume fraction,as well asthe local Machnumber.

The currentresultsprovideevidence thatat10% volume fraction,

Fig. 7. Impulse for all particles over the full simulation time t end≈113 τp. Circles:

particle-resolved data. Triangles: Spatially averaged particle-resolved data, black squares: spatially averaged isolated particle model, light gray squares: spatially averaged velocity corrected model, dark grey squares: spatially averaged volume fraction corrected model.

theshock-diffraction patternsaresuﬃciently disruptedinaman- ner thatremovesthe negativeforceon average.However, further studies are warranted in order to characterize how thisdepends onbulkﬂowproperties.

For prediction of particle movement, the impulse is a better indicator than the force imposed at any speciﬁc time. Thus the impulseisa goodquantityto considerwhenevaluatingthe overall performanceof thedraglaws.Withtheparticle-resolved simulation data as input, the volume-fraction corrected model gives the best prediction of impulse, except in the downstream edge region, where the velocity corrected model is better. In Fig. 7, the impulse over the whole simulation time is shown for each particle, along with the impulses predicted by the drag models.

The particle-resolved data is shown for each individual particle as well as spatial average over bins with width L/60. The impulses predicted by the drag models are only shown as spatial averages. The isolated particle model underpredicts the average impulse by up to 50%. The comparisons with the velocity correctedmodel and thevolume fractioncorrected models are better,buttheystillunderpredicttheimpulse.Thevelocity-corrected model isbetter at thedownstream edge due tothe effect ofin- creasedMach number causedby increased velocity.As shownin Osnesetal.(2019a),theincreasedforcesatthedownstreamedge in the particle-resolved simulations are due to a rapid increase in the Mach number. The drag forces depend strongly on Mach numberin thetransonic regime, see e.g. Nagataet al., 2016,and thedraglawofParmaretal.(2010)accountsforthisdependence.

Thus,thevelocity correctionismoreimpactfularound thedown- streamedgeduetothehigherlocalMachnumbershere.Thedrag law of Parmar et al. (2010) captures the increased forces at the downstreamedge,andthishasbeenshowntobeparticularlyim- portant for simulations of shock-particle cloud interaction, since accumulation,ratherthan dispersion,of particlesnearthe down- streamedgehasbeenshowntobeaconsistentprobleminbothEE andELsimulations(TheofanousandChang,2017).Still,theconsis- tentunderpredictionofimpulseevenwithproperlyMachnumber dependentdraglawsindicatesthatsimulationsthatusethesedrag lawsare likelytounderpredict particlemovement anddispersion duringtheshock-particlecloudinteraction.

The unsteady forces,andthe corresponding impulses, are im- portantonlyforashorttimeaftertheshockwavepassesovereach

(9)

Fig. 8. Impulse as a function of time for the different force components. (a): Isolated particle model. (b): velocity-corrected model. (c): volume-fraction corrected model.

particle.Fig.8showstheimpulse J=

t 0

p (depending onthe model).As time goeson, the impulse due to both unsteady forces becomes lessandlessimportant.The ﬁgureshowsthat onlyparticles that have response time-scales of less than O

10

τ

^p

will havea mo- tionthatisnoticeablyaffectedbytheinviscidforces.However,the responseofthegastotheseforcescanstillbesuchthattheinclu- sion ofunsteadyforces isimportanteven forparticleswithlarge responsetime-scales.

There are very small differencesin the relativeimportance of the force terms for the three models. The most noticeable dif- ference is that the pressure gradient impulse is more signiﬁcant for the isolated particle-model when compared to the velocity- corrected and the volume-fraction corrected models at the later stages.

3. Comparisonofparticle-resolvedandELsimulations 3.1. ELSimulationapproach

ThegoverningequationsfortheELsimulationsarethevolume- averaged mass, momentum and energy conservation equations.

Since the problem under consideration is one-dimensional, only the corresponding one-dimensional volume-averaged equations willbegivenhere.Additionally,theparticleswillbeassumediden- tical, stationary,andinert. Inthefollowing,

ψ

^denotes^a^phase-

averaged quantity, where

ψ

îs ânyfluid quantity, and is related to the volume-averaged value,

ψ

, through

α ψ

=

ψ

^. ^Addition-

ally,

ψ

˜=

ρψ

/

ρ

^denotes ^a Favre-averaged quantity. The devi-

ationfromtheFavre-averaged quantitywillbe denoted

ψ

)

−

∂

^x

α ρ

ûûû^˜i

+

ξ

, (25)

where

μ

^is^the^dynamic^viscosity, ^E^˜=e˜+¹₂u˜²+¹₂uu isthe to- talenergyper unit mass,whereeis theinternal energyperunit mass,and

λ

^is^the ^heatconductivity,which isassumedtobe related to the viscosity through a constant Prandtl number of 0.7.

NotethatthestationaryparticlesdonoworkintheEulerianframe, and thus there are no drag-related terms in the energy conser- vationequation. Intheenergyconservationequation,mostofthe sub-gridscaletermshavebeencollectedintheterm

ξ

,whichwill besetto0inthiswork.Theimportanceofthesetermsisnotwell known for shock-wave particle cloud interaction, and to the au- thorsknowledge,noappropriate modelsexistforthesetermsun- derthecurrentflowconditions.Theomittedsub-grid scaleterms can be foundin Osnes(2019). Mehta etal.(2020) quantified the inviscidsub-gridtermsusinginviscidparticleresolvedsimulations, and found that both the internal energy-velocity correlation and the pressure-velocity correlation were negligible. The remaining sub-gridtermsarethethirdvelocitymomentandtermsinvolving themass-weightedturbulentvelocity,whichhavepreviouslybeen modeled with an additionaltransport equation Schwarzkopf and Horwitz (2015). The mass-weighted third velocity moment can be neglected based on the principle of receding influence. In VartdalandOsnes(2018),theproductionoffluctuationkineticen- ergy bymass-weighted turbulent velocity termswasfound tobe smallincomparisontotheotherproductionterms.Basedonthis, itcan beassumedthat themass-weighted turbulent velocitydif- fusiontermsinthetotalenergyequation arealsosmall,andthus

(10)

theywillbeneglectedhere.Anin-depthassessmentoftheimpor- tanceofthesetermsisatopicforfuturestudies.

Onlythevolume-fractioncorrectedmodelwillhaveanon-zero valueforuu,whichisobtainedbyassumingthatuu≈

^u^u

andusingEq.(16).Thismeans thatthe densityﬂuctuation isnot correlated withthe magnitude of the velocity ﬂuctuation, which is unlikely in general, buta moredetailed modelis necessaryto accountforthiscorrelation.

The equation of state forthe gas is the ideal gas law, where internalenergy,pressureanddensityarerelatedby

p=

( γ

⁻¹

) ρ

^e, ⁽²⁶⁾

with

γ

=1.4. The temperature, T is related to the internal en- ergybyaconstantheatcapacity,andweassumethattheviscosity varieswithtemperatureas

μ (

^T

)

⁼

μ

ref

_T

T_ref

⁰.76

, (27)

where

μ

ref isthevalue oftheviscosityatthereferencetempera- tureT_ref.

The conservationequationsEqs.(23)to (25)aresolvedwitha control-volume based finite volume method. The fluxes between control volumes are computed with a modified HLLC Riemann solver,andMUSCLreconstructionwiththeminmodlimiterisused todefinetheRiemannproblems.Thesolutionisadvancedintime withathird-orderRunge-Kuttascheme.Thefluidvariablesarein- terpolated linearlyto the particle positions. Severalinterpolation pointsareusedforeachparticle.Thesepointsareassignedweights depending on their position relative to the particle center, and the sum of the interpolated valuesto each of thesepoints mul- tiplied bythepoint’s weightapproximates thevolume-average of the fluid property over the particle volume.The force fromeach particle on thefluid is transferred witha filterto ensureconsis- tency upon mesh refinement. We use the two-step filtering ap- proachofCapecelatroandDesjardins(2013),whereamollification isfirst usedtotransfertheLagrangian datatotheEulerianmesh, followed by a diffusion step to obtain the desired spatial distribution ofthetransferredterms.Thefilterwidthrecommendedin CapecelatroandDesjardins,2013was

σ

f=3Dp,whichiswhatwe usehere.

We use a four times finer spatial resolution in the EL simulations than the sampling-resolution of the volume-averaged particle-resolveddata.ThespatialresolutionfortheELsimulations isthereforeapproximately^x=0.125Dp.Thegridresolutiondoes haveaninfluenceonthesimulationresults,andgrid-convergedre- sultscan notbe expectedatleastuntilthefilterkernels arewell resolved. The dependence of the simulation results to grid reso- lutionis showninAppendixA.Sincethe objectiveofthecurrent work isto evaluate the applicability of the drag-laws, we accept spatial resolution that is likely not achievable in simulations of manyreal-worldproblems.Thesesimulationswilloftenbeunable toresolve anyflowattheparticlescale.However,onlyafewfea- turesintheresultsareverysensitivetothegridresolution.These arethespeedatwhichthereflectedshockwaveisgenerated,and theflowexpansionatthedownstreamparticlecloudedge.

3.2. Simulationresults

In general, the mean flow fields are quite well captured with all threeapproaches, butthe velocity corrected model is slightly better thanthe others.Fig. 9showsthemassdensityatfourdif- ferenttimesfortheELsimulationsandtheparticleresolvedsim- ulations.Earlyon,theELsimulationsagreewellwiththeparticle- resolved simulations. The strength ofthe reflected shockwave is well predictedwith the velocity-corrected modelearly on,while the twoother modelshavea slightlyweaker reflectedshock.The

volume-fraction corrected modelgives a stronger reflected shock thanthe isolated particlemodel,asexpecteddueto thestronger forces imposed on the flow by the particles. The stronger shock reflectionwith the velocity-correctedmodel can be attributed to the velocity-fluctuation term, which imposes an extra upstream- directed force on the flow at the upstream particle cloud edge.

Similarly,anadditionalforceappearsatthedownstreamedge,this timedirecteddownstream.Thisleadstoastrongerﬂowexpansion atthedownstreamedge,whichagreeswellwiththeoneobserved in the particle-resolved simulations. The isolated-particle model andthevolume-fractioncorrectedmodelbarelypredict anyover- expansionatall,andthustheagreementinthedownstreamregion becomespooratthelatertimepoints.Itshouldbenotedthatthe expansionregionisparticularlysensitivetothegridsize,asshown inAppendixA.

The corresponding pressure and velocity ﬁelds are shown in Figs.10and11.Thevelocityresultsofthevelocitycorrectedmodel isslightlyhigherthantheparticle-resolvedvelocityﬁeldintheex- pansionregion att=60

τ

^p^. În Ôsnes êt âl.^(2019, ²⁰²⁰⁾^, ^the âu-

thorsfoundthatthevelocityfluctuationsdecreasedrapidlyatthe downstreamedge,butnottozero,andnotassharplyaspredicted by the velocity-corrected model. The too sharp velocity fluctuation decaypredictedby themodel produces astrong streamwise force that appears immediately after the shock wave passes. In the particle-resolvedsimulations, there is a slighttime-delay be- foresignificantvelocityfluctuationsappear,anditisthereforenot surprising that the flow expansion is stronger early on with the velocitycorrectedmodel.

Fig.12showstheaverageparticleforcesinthe ELsimulations alongwiththeaveragedforcesintheparticle-resolvedsimulations.

Fig. 13shows thecontribution ofthe different force components intheELsimulations.Inthesesimulations,likeinthecasewhere weusedvolume-averaged particle-resolveddataasinput, weﬁnd that the peak force is overpredicted by the volume fractioncorrectedmodel. Here, thisis alsothe case forthe isolated particle model,butthe velocity correctedmodel underpredictsit instead.

Thetimeofthepeakforceisalsoslightlydelayedcomparedtothe one intheparticle-resolvedsimulations. However, sincetheforce historyforeachparticleis shiftedwitht0,p,thedefinitionofthis time-pointhasadirect impacton thelocation ofthepeak value, especially when the computational grid is not very fine, so that thenumericalsmoothingoftheflowfieldsleadstoanearliert0,p

fortheELsimulations.Nevertheless,theshockrelatedtransientis signiﬁcantly smoothedin time,andthe forcesdo notincrease as quicklyfollowingshockimpactasfortheparticle-resolvedsimulations.FinerEuleriangridswillhaveaneffectonthis,andwilllikely lead tomorerapid force increases.However, thisisalsolikely to increase the peak value further, and thus bring the peak forces further away from those observed in the particle-resolved simulations.Thevolume-fractioncorrectedforceisstillthehighest,but thedifferenceissmallerthanitwasintheresultspresentedinthe previoussection.

At late time, the forces with the volume-fraction corrected modelareagainthosethatareclosesttotheparticle-resolvedsimulations, but they are in less agreement than when the forces were evaluated with the volume-averaged particle-resolved flow fields. This is related to the differences in the mean flow fields, c.f. Figs. 9, 10, 11. The velocity-corrected model resultsare simi- laratlatetimesastheywere withthe volume-averagedparticle- resolved data as input, while the results of the isolated particle modelaresurprisinglyslightlybetterthantheywerepreviously.

Considering the force components,the largest differences between thethreemodels are foundin thequasi-steadyforcesand theundisturbedﬂowforce.Thequasi-steadyforceislargestforthe velocity-corrected model early on, butis largest for the volume- fractioncorrected model atlate time. It issigniﬁcantly lower for

(11)

Fig. 9. Density proﬁles at four times. (a): t = 15 τp, (b): t = 30 τp, (c): t = 60 τp, (d): t = 100 τp.

Fig. 10. Pressure proﬁles at four times. (a): t = 15 τp, (b): t = 30 τp, (c): t = 60 τp, (d): t = 100 τp.

theisolatedparticlemodelatalltimes.Theundisturbedﬂowforce islowerearlyonforthevelocity-correctedmodelthantheisolated particle model and thevolume-fraction corrected model. Forthe isolatedparticlemodel,thiscanbeexplainedbytheweakershock waveattenuationwhichleadstoastrongershockwaveandthere- fore stronger pressure forces. For the volume-fraction corrected

model, the increased undisturbed ﬂuid force is more surprising, sinceitistheonlyforcecomponentthatisnotincreasedbyascal- ing factor. It is a resultof an increasedpressure gradient, which isinducedbylargerforcesandnocorrectionduetovelocity ﬂuc- tuations in the equation of state. The inviscid unsteady force is

(12)

Fig. 11. Velocity proﬁles at four times. (a): t = 15 τp, (b): t = 30 τp, (c): t = 60 τp, (d): t = 100 τp.

Fig. 12. Average particle force predictions for all particles in the particle-resolved simulations (black lines), and for the EL simulations (colored lines). The shaded areas indicate one standard deviation. (a): Isolated particle model, (b): velocity corrected model, (c): volume fraction corrected model.

alsoslightlyhigherinthevolume-fractioncorrectedmodel,which againcanbeattributeddirectlytothescalingfactor.

Fig.14showstheimpulseoverthefullsimulationtimeforthe particleresolveddataandthethreeELsimulations.Thefigurealso showshowthetotalimpulse(sumoverallparticles)developsover time. Consideringfirst how the impulse atlate time varies with position,wefindasexpectedthattheagreementwiththeparticle- resolved data is worst for the isolated particle model,while the velocity-corrected and the volume fraction corrected models are

thebestinthe edgeregions andthecentralregions, respectively.

It isworth notingthat all three models predicta larger impulse, and thus higher velocity, for the particles near the downstream edgethan thosefurtherin.Itwasemphasizedin Theofanousand Chang(2017)thatthischaracteristiciscrucialtocaptureinsimu- lationsofshock-accelerated particlelayers.TheMach-numberde- pendent quasi-steady dragmodel ensures that this characteristic is indeed captured by the EL simulations. The impulse is however about50% higher for the velocity-corrected model than for