Mesoscale simulations of shock compaction of a granular ceramic: effects of mesostructure and mixed-cell strength treatment

The shock response of granular materials is important in a variety of contexts but the precise dynamics of grains during compaction is poorly understood. Here we use 2D mesoscale numerical simulations of the shock compaction of granular tungsten carbide to investigate the effect of internal structure within the particle bed and ‘stiction’ between grains on the shock response. An increase in the average number of contacts with other particles, per particle, tends to shift the Hugoniot to higher shock velocities, lower particle velocities and lower densities. This shift is sensitive to inter-particle shear resistance. Eulerian shock physics codes approximate friction between, and interlocking of, grains with their treatment of mixed cell strength (stiction) and here we show that this has a significant effect on the shock response. When studying the compaction of particle beds it is not common to quantify the pre-compaction internal structure, yet our results suggest that such differences should be taken into account, either by using identical beds or by averaging results over multiple experiments.


Introduction
Despite the ubiquity of granular materials on Earth and throughout the solar system their response to shock is poorly understood, particularly their internal dynamics. These processes are vital for understanding large-scale shock events, such as asteroid collisions and detonations, and are instrumental in defending against them. A great deal of experimental research into the bulk response of granular materials, especially ceramics and sand, currently exists [1][2][3][4][5][6]. However, the high strain-rates and short timescales involved in shock compaction have made experimental work into granular dynamics challenging. Consequently, explicit investigation of grain dynamics has tended to employ 'mesoscale' computational simulations using Eulerian shock physics codes, such as CTH or ALEGRA [7][8][9]. Such simulations explicitly resolve hundreds-to-thousands of grains, such that microscale effects are present but the bulk response is also considered. This direct link between the micro-and macroscales has proven instrumental in furthering understanding of granular dynamics in shock compaction [10][11][12][13][14].
Several recent mesoscale simulation studies have investigated the effects of specific microscale phenomena, such as sensitivity to intergranular shear resistance [11,15], grain morphology [16] and grain fracture [17], on the macroscale shock response of particle systems. However, as more microscale effects are implemented there is also a need to consider the effects of certain model assumptions that are poorly tested [11,18]. In particular, it is common in mesoscale simulations to use either a single grain arrangement or different random arrangements when constructing a synthetic Hugoniot curve for the particle system. Yet, grain arrangement has been shown to have a discernible effect on shock response [11,[19][20][21] with some studies directly investigating the effects of varying bed structure [18,21]. In general, these investigations quantify variation in structure between particle beds by the exposure to their respective bed-generation algorithms. They make no attempt to analyse the final bed structure and can not guarantee that their beds do not exhibit structurally similar features.
By way of contrast the effect of grain structure on quasi-steady state deformation is well studied and known to be significant [20,[22][23][24][25]. A number of metrics, such as the coordination number [26], have been developed to quantify and track the mesostructure in a particle bed. The coordination number is generally defined as the average number of contacts that an entity has with other entities in the system [24,26]. These entities are typically particles but might instead represent voids [26]. In this work, which considers two-dimensional beds of circular particles, the coordination number (which we label as 'Z' in accordance with its use in the literature [24]) represents the average number of independent points of contact on the circumference of each disc.
Here, we examine the effect of coordination number Z on the planar shock compaction of mono-disperse, 2D, particle beds using mesoscale numerical simulations. Our simulations are based on previous laboratory and numerical experiments of the shock response of granular tungsten carbide [27][28][29] to provide a means of benchmarking iSALE simulation results with previous work. The porosity (and hence initial bulk density) of a particle bed will have a significant effect on both the shock response and the range of possible coordination numbers.
We therefore kept porosity constant between comparative simulations so that particle beds differed only in Z. We show the mesostructure within the particle bed directly affects the shock response of the particle bed: an increase in coordination number ('more' structure) shifts the Hugoniot (locus of all shock states) to higher shock velocities, lower particle velocities and lower shocked densities for a given impact velocity. In addition, a sensitivity to structure implies a sensitivity to shear resistance between grains (approximated using stiction) and we also show that removing intergranular resistance (no stiction) drastically reduces the resistance to compaction. Most simulations used the iSALE shock physics code [30] and several additional simulations were performed with CTH [7] for comparison. The good agreement between the codes demonstrates that the results are robust within Eulerian frameworks and not specific to either code.

Setup and methods
Numerical mesoscale simulations were designed to simulate a 1D flyer plate impact experiment. They were loosely based on previous 2D planar mesoscale simulations performed by Borg and Vogler [27] that modelled the compaction of a granular tungsten carbide (WC) bed with a porosity of 45% [1]. This study considers an aluminium (Al) impactor, with an initial velocity (V I ), that strikes a particle bed of tungsten carbide (WC) with 50% porosity (f, defined as the ratio of void space to occupied space in the particle bed). No driver plate was used. The full setup is presented in figure 1. The impactor measured 21×1 mm, however only the millimetre of the impactor closest to the particle bed was resolved at high resolution; the remainder of the impactor was modelled with progressively larger cell sizes away from the particle bed, to minimise the computational cost of modelling a long impactor. A long impactor was necessary to delay the arrival of the release wave from its rear until after the compaction wave had traversed the particle bed. The particle bed measured 1×1 mm. A range of resolutions were considered, from 10 to 20 cells per particle diameter (cppd), and these were offset against their relative computational costs. Consequently, a resolution of 2 μm/cell was used and a particle resolution of 16 cppd employed. This gave the particles a characteristic diameter of 32 μm. Previous studies found that particle size and morphology variations did not affect the shock response of the simulated particle beds [18]. Hence, a mono-disperse size variation was employed here to simplify and constrain the scope of the study, as was done in previous work [27]. The boundaries of the simulation, perpendicular to impact, were frictionless, rigid and not periodic. Open, 'outflow', boundary conditions were used at the top and bottom (behind and ahead of the impact, respectively). The aim of these simulations was to quantify the shock response of the particle bed up to the pressure required for complete compaction. The shock response is typically represented as a 'Hugoniot' curve: this is the locus of all possible shock states for a material. To do this, a series of simulations were performed using a range of impactor velocities from 150-950 ms −1 , in intervals of 50 or 100 ms −1 , resulting in Hugoniot curves with approximately 10 points.
A Mie-Grüneisen equation of state (EoS) was used for both the aluminium (Al) and the WC as was done in previous studies [18,27]. The associated constants are presented in table 1. A Johnson-Cook strength model was used for the Al [31], while the WC particles were treated as elastic, perfectly plastic with a constant yield strength, independent of any other variables such as strain and temperature. Neglecting thermal softening of the WC (and use of a simple Mie-Grüneisen EoS) is justified given WC's high melting point (3143 K). Even at the highest impact velocity considered the bulk shock temperature of the WC was ≈1200 K, with a standard deviation of ≈400 K. Only a small percentage of cells within the  Figure showing the geometries used. Simulations were performed in the horizontal plane; we refer to the x and y directions as longitudinal (L) and transverse (T), respectively. The mesh is Cartesian, and structured and the simulations were not periodic in the T direction but used solid 'freeslip' boundary conditions. The extension zone used a decreasing resolution to minimise the computational cost of simulating a large impactor. The porosity is denoted by f. The iSALE shock physics code [30] was used to perform most of the simulations in this work and has been used in the past to simulate mesoscale compaction of primitive solar system materials [32][33][34]. iSALE is an Eulerian code with a fixed mesh within which particles move. When particles come into contact their shared interface is represented by so-called mixed cells that are a mixture of each particle material (and void). Whilst iSALE does track the interfaces between materials in a mixed cell it does not track the interfaces between particles. Hence, the intergranular resistance to shear is represented by the strength of these mixed cells: an approach referred to as 'stiction' [11]. iSALE can distinguish between a number of different materials using volume fraction fields, regardless of whether the properties of these materials differ. We refer to these different material volume fraction fields as material numbers. In the simulations five different material numbers were used to represent particles in the bed, with assignments designed to minimise the number of contacts between grains of the same material number. Such instances result in unrealistic artificial 'welding' of particles. All particle beds that were used had zero contacts between particles of the same material number at time zero; five material numbers was the minimum required to achieve this. It is possible for the number of contacts a particle has to exceed four during a simulation; however, at such a point the distinction between particles becomes less important because of the removal of most of the porosity.
This study investigates the effect of two parameters on the granular WC Hugoniot: the initial arrangement of grains in the particle bed, and the treatment of strength in cells containing multiple materials (grain-grain contacts). Mesostructure was quantified using the coordination number, Z, the average number of contacts each particle has with other particles. The initial density (and hence, porosity f) was kept constant between beds comparing the effects of changing Z. In a particle bed with a high Z, the network of contacts typically covers more of the bed, and should provide a greater resistance to deformation. However, mesostructure relies on resistance between grains to have an effect. This resistance is approximated in (Eulerian) shock physics codes by stiction: the strength assigned to cells containing multiple materials. Two mixed-cell strength treatments were compared in simulations: in the first the strength of a cell was calculated as a weighted average of the individual strengths of each material occupying it (finite stiction); in the second (the 'friction-less' case) cells with multiple materials were automatically given a strength of zero (no stiction).
2.1. Quantifying mesostructure in the particle bed Multiple particle beds, containing randomised arrangements, were generated for the purposes of these simulations. Their bulk porosity f and approximate number of particles were kept constant (to within ≈ 0.01%), along with the particle size and circular shape. Gravity was not included in simulations to allow for the generation of more varied values of the coordination number; that is, particles could be 'suspended' with zero contacts. Whilst this is an unlikely occurrence in a physical ensemble under terrestrial gravity, such arrangements can occur in very low-gravity environments. Such an arrangement also approximates the distribution of particles in a 2D slice through a 3D particle bed or in beds where the gravity vector is perpendicular to the impact direction.
The amount of mesostructure in the bed was controlled by varying the coordination number Z. The average number of contact points in the particle bed is heavily dependent on the initial porosity of the particle bed, the definition of a 'contact', and the algorithm used to place particles. Any one algorithm can have a multitude of potential particle arrangements, each with significant variation. A bed created at a given initial porosity with a fixed number of particles can only contain a finite number of possible Z values. A theoretical maximum value of Z = 6 is possible for a 2D bed of mono-disperse discs at the maximum packing density (MPD, ≈10% porosity, ≈90% volume fraction). This is for circles in an infinite plane; in a finite rectangle the true value will be slightly more porous and is equivalent to the theoretical maximum density for that domain. The effect of mesostructure will be more visible in particle beds with initial densities significantly below maximum packing, which are capable of producing the full range coordination numbers, Z 0 6  < .
Each particle was considered individually to calculate its number of contacts (figure 2) and the mean was calculated. In the computational domain, each particle is represented by a circular region of cells. At time zero there are no mixed cells; all cells inside a particle are full. A 'contact' between two particles at time zero in the simulation was defined as two, or more, adjacent cells (diagonals were not considered adjacent) from two different grains with no void cells separating them. For the particle beds in this work, we generated a library of beds with constant porosity f and a range of Z by varying the particle placement algorithm's input parameters and selecting those beds that fell within a pre-defined range of Z (figure 3). The algorithm employed here placed particles into the bed in two ways. A proportion of the total number of particles (one of the input parameters which was chosen beforehand) were inserted at random into an empty space in the domain, not in contact with any other particles and immediately fixed. The remainder of the particles were first inserted in the same manner as previously described, but were then moved in a random direction until they made contact with another grain. Varying the ratio of grains that were placed versus placed-and-moved allowed for variations in Z. Only grain centres were required to be within the boundaries of the bed and so partial particles were allowed on all boundaries.
Desired bed arrangements were those with variable Z but constant and approximately spatially uniform porosity, f. However, for high coordination numbers some beds displayed significant porosity gradients (figure 3 shows an exaggerated example of such a bed that was Figure 2. A subsection of a particle bed; contacts between particles are denoted by white lines, and the white numbers denote the number of contacts for that particle. Only non-zero values are shown. not used). To distinguish (and reject) such beds a second metric was used: the maximum porosity difference across the bed, K. This was calculated by first dividing the bed into five equal columns and calculating the porosity in each column. The maximum difference between the columns was recorded. This calculation was repeated along the orthogonal axis of the bed (rows in this case) and K was defined as the largest difference in porosity between any two rows or columns. Among the 272 particle beds generated, the maximum porosity difference across the bed was found to be small and relatively constant K 0.1  for beds with a coordination number Z 1.9  (figure 3). Hence, only particle beds with Z>1.9 were excluded from the study. Figure 4 shows two examples of particle beds that were chosen for this study. None of the generated beds exceed Z = 2.5 or come close to the maximum value Z=6, because of the random nature of the particle placement algorithm and the interdependence of f, Z and K. There is negligible chance of constructing a bed with Z≈6 at random and with a bulk porosity of 50% such a bed would necessarily have large K and be excluded.

Calculating the Hugoniot
To calculate the state of the particle bed behind the shock for a given simulation it was necessary to extract a value of shock velocity (U s ), particle velocity (u p ), bulk density (ρ) and bulk longitudinal stress (σ LL ). The shock velocity is the speed of the shock wave's movement through the particle bed, which is calculated as the slope of a linear regression of shock front position against time. The other variables are calculated as averages over the particle bed behind the shock front of longitudinal material velocity, mass density and longitudinal stress. In this work, we used two different methods to identify the shock front position and calculate the shock variables for constructing the Hugoniot (the 'transverse average' method and the 'stress contour' method, figure 5, right and left respectively).
The first method produced a 1D wave profile in the longitudinal direction at each time step by transversely averaging all fields across the granular bed (figure 5). We refer to this approach as the transverse average method. The 1D profile was used to locate the instantaneous shock front position and find the shock state behind the shock front. The shock front position was identified as the longitudinal position at which the magnitude of the longitudinal stress ( S P LL LL s = -, where the deviatoric stress, S LL , is negative in compression, and pressure, P, is positive) exceeded half its value behind the shock front from the previous time step. (At time zero the longitudinal shock stress was taken to be zero.) The instantaneous shock state (i.e., particle velocity, density, stress) at each time step was found by calculating the mean value of the relevant field between the shock front position and the impactor/particle bed interface position. Cells with no mass (void) were included to calculate the average shock stress and mass density, but not particle velocity. The time-averaged shock state was found by averaging the instantaneous shock states over all times whilst the shock remained 'steady'. This period was determined to begin after the shock had travelled two particle diameters and lasted until it reached 0.85 mm along the bed. This second threshold ensured that no reflections from the end of the simulation contributed to a shock state. This transverse averaging method has been used in multiple previous studies, including those using CTH [18,27]. A second method to locate the shock front used a contour of longitudinal stress (figure 5) [34]. We refer to this approach as the stress contour method. The selected stress contour was chosen to be equal to half the value of the instantaneous longitudinal shock stress value from the previous time step. The contour is undefined at time zero. To track the shock front position with time, we defined the shock position as the median longitudinal position of the shock front contour. The instantaneous shock states of density, longitudinal stress, and particle velocity were found by taking the mean of their values between the shock front contour and the interface between the particle bed and the impactor. Again, empty cells were excluded when averaging particle velocity. As with the transverse average method, the time-averaged shock states were found by averaging the instantaneous shock states over the period that the shock wave was steady. The difference between the transverse average and stress contour methods is illustrated and discussed in section 2.5.

Benchmarking iSALE and CTH
Before using iSALE to explore the effect of mesostructure on the shock compaction of particle beds, several simulations were run with both iSALE and the CTH shock physics code [7] to examine intercode variability. Simulations to benchmark the two codes were set up as close as possible to each other, which required some minor modifications to the setup described above. The benchmark simulations used the same material properties, geometries, boundary conditions (non-periodic), mixed cell strength treatments and the same particle arrangement ( figure 6). However, for this comparison a bulk porosity of f = 45% (≡ 55% volume fraction) for the particle bed was used to allow an additional direct comparison with the experimental data [27]. The benchmark simulations also used slightly smaller grains with 12 cppd, and diameters of 24 μm because this was the preferred resolution of CTH. The transverse average method was used to calculate the Hugoniot states from both the iSALE and the CTH simulation results. The impactor position was chosen as the position of the impactor/particle bed interface that was furthest along, longitudinally.

Results of the iSALE-CTH comparison
The Hugoniot response (in U s -u p space) inferred from the CTH and iSALE simulations using identical particle arrangements and mixed-cell strength treatments are in good agreement (figure 7). Differences in shock and particle velocities between the simulations is consistently 4%  ( figure 7). Previous work has suggested that the granular tungsten carbide Hugoniot consists of approximately four linear 'sections' (in U s -u p space), possibly reflecting different regimes of behaviour [18,27,29]. Although the data presented here do not extend to   [27]. (Right) Percentage difference between the iSALE & CTH shock and particle velocity results. Only comparable data is presented here. The two results never differ by more than 4% for either variable.  Figure 7 also includes published experimental Hugoniot data for a WC powder [27]. The characteristic grain size (32 μm) of the WC powder used in the experiments was larger than that used in our simulations and the grains were highly angular. In other regards, the experimental setup was very similar to the simulation setup used here. Over the limited range in particle velocity for which experimental data exist, the experimental and numerical data show a very consistent U s -u p slope (experiment≈1.46; iSALE≈1.56; CTH≈1.32), but the numerical data lie somewhat below the experimental data. We attribute this offset to the 2D approximation used and a potential underestimation of the strength of individual WC particles. The bulk strength of the particle bed may be underestimated because of the lack of out-of-plane rigidity produced by a 3D grain network. Previous studies have found that fully 3D simulations, with stiction, over-predict the shock response, shifting the Hugoniot too far up in U s -u p space [18]. 3D simulations without stiction, however, under-predict the shock response; they were found to be in good agreement with 2D simulations with stiction [18]. These findings suggest that the use of 3D would address the under-prediction of the shock response, however, they would not be sufficient to bring the results in line with experiment. Alternatively, the simulations may not accurately represent the strength of the WC particles. The dynamic yield strength of WC is not well understood and is likely strain and strain rate dependent. In this work, the shear strength of WC is treated as a constant (see table 1), with an estimated value derived from a measurement of the Hugoniot elastic limit of 6.6 GPa [35]. Moreover, previous work using CTH found that increasing the yield strength of the individual grains from 5 to 7.5 GPa was sufficient to bring the numerical results in line with the experimental Hugoniot [18,27,29]. Hence, it is possible that uncertainty in the WC yield strength or an oversimplification of the WC strength behaviour account for the disparity between experiment and simulation.

Comparison of shock state extraction methods
To compare and contrast the two methods for extracting shock state variables, described in section 2.2, from our mesoscale simulation results, we calculated the instantaneous and timeaveraged shock states for a number of simulations using both methods. The two methods differ in their approach to identify the shock front and in how they calculate instantaneous shock states behind the shock front. While the two methods produce similar results, differences in shock velocity of up to ≈10% (figure 8) and particle velocity of up to ≈15% were observed using the two approaches. In many cases, a steady shock wave was identified much earlier in the simulation using the stress contour method compared with the transverse average method (figure 8). For this reason, the stress contour method produced more consistent timeaveraged shock states and Hugoniots with significantly less 'scatter'.
The efficacy of the transverse average method depends on the 1D wave profile of the shock being an accurate representation of the actual shock, which is compromised when the shock wave is non-planar. In such cases, the 1D approximation will exclude some shocked material from the shock state calculation, and similarly, will include some unshocked material ( figure 5). The closer the shock front is to a planar wave, the smaller this effect will be and the better the transverse average method becomes. The second, stress contour method does not have this limitation because it accounts for the 2D variation in shock front position. It is able to better represent the heterogeneity of the shock front that is characteristic of shock propagation in particulate systems at the mesoscale. For these reasons, the stress contour method was used to calculate the shock states in the mesostructure investigation.

Sensitivity of the Hugoniot to mesostructure
Hugoniot curves in U s -u p space were constructed for suites of simulations with the same initial grain arrangement and coordination number Z ( figure 9). As all post-shock parameters can be affected by changing Z, only the initial conditions (V I , Z) are considered to be independent variables. A common trend in the results is for the U s -u p slope to be steep at low u p , shallower at intermediate u p and then steepen again at higher u p . There is no sharp discontinuity between these regions and the Hugoniots appear to transition smoothly as u p increases. At high particle velocities (>400 ms −1 ) we observe the clearest variation with Z, as u p decreases the existence of a monotonic relationship between U s or u p and Z becomes more unclear.
Subtle differences in the Hugoniot curves are difficult to discern in U s -u p space because it is not particularly sensitive to different material properties [36]. Shock velocity in the moving frame, U F ( U u s p = -), as a function of u p , is generally more sensitive to subtle variations in material properties [36]. In this space, the differences and variations observed in U s -u p space are amplified and the change in Hugoniot shape with respect to Z is more pronounced (figure 9). At particle velocities exceeding 400 ms −1 there is a systematic shift in the Hugoniot to lower u p and higher U s with increasing Z at constant V I . Values of σ LL and ρ were extracted explicitly, as described in section 2.2, rather than using the Rankine-Hugoniot jump conditions. However, unlike u p , these fields exhibited notable variation across the particle bed behind the shock (figure 10). In particular higher stresses were observed along force chains in the bed whilst grains that did not take part remained largely unaffected. The level of variation we observed is consistent with that presented in the previous studies upon which this work is based [27].
The Hugoniots in stress-density space have less variation in the y-axis (stress) than in U u s p -space, with all points achieving approximately the same stress for each impact velocity Figure 8. Shock position versus time (left) and particle velocity u p versus shock position (right) for the 'Transverse Average' and 'Stress Contour' methods. Shock positions more than six particle diameters into the bed and less than 0.85 mm (solid horizontal and vertical lines respectively) from the end were used to calculate bulk values of shock and particle velocity (dashed pink and black lines, respectively).
( figure 11) although there appears to be a slight negative correlation of stress with Z for high V I . In this space density is more sensitive to coordination number (mesostructure), and the Hugoniot shifts to lower densities as Z increases. However, the difference decreases for higher stresses to the point it appears to vanish as they intersect with the fully dense WC Hugoniot. Above this point the particle bed will achieve full compaction and removal of all void space regardless of the initial mesostructure, and the bulk density will tend to approximately the solid density (≈15.56 g cm −3 ). Figure 12 compares two suites of iSALE simulations that use the same initial particle bed (Z = 1.5) but differ in their treatment of strength in mixed cells; 'with stiction', a weighted average of the strengths of material in a cell (which was used for all other simulations so far presented in this work) and 'without stiction', i.e. 'frictionless' (any mixed cell has zero strength). This will make almost all structures involving particles 'unstable' and liable to break down until the particles reach their most stable arrangement possible in the time allowed.

Sensitivity of the Hugoniot to the strength treatment of cells containing multiple materials
The weighted average results appear as before (figures 9 and 11), but removing stiction shows a single, approximately linear, trend in U F -u p space. The slope of this linear trend is steep throughout and is most similar to the slope of the high-u p section of the Hugoniot for the weighted average case, where the two intersect with each other ( figure 12). The removal of stiction allows particles to 'slip' past each other which removes the effects of mesostructure in the particle bed. In LL s | |-ρ space the particle bed is compacted to densities exceeding the MPD (for square packing f≈22%) for all data points because particles are free to move until they reach a stable position. Similarly the stiction-less Hugoniot intersects with the weighted average Hugoniot, and the Hugoniot for initially solid WC, for high ρ. -)-u p space. This space amplifies differences between Hugoniots and presents a more accurate representation of the shock response [36]. Standard errors were calculated for each Hugoniot point, however they were not included because they were too small to be visible.
The intersection between these two Hugoniots, and the solid WC Hugoniot, indicates that in this regime the shock response is insensitive to stiction between grains. This is as expected; when the V I is large enough to compact all porosity out of the bed the effects of grain-grain interactions on removal of void will no longer matter. Above this point there is no remaining porosity and the response of the homogeneous material dictates how the system evolves. Additionally, this result implies that despite the presence of porosity, the temperature of the shocked particulate matter is not sufficient to have a noticeable effect on its response. If this were not the case, the porous Hugoniot would not intersect with the 'cold', solid, WC Hugoniot.  . The shock response of the particle beds in figure 9 in stress-density space. Dark to light colours represent increasing values of Z. The theoretical 'maximum packing density' (MPD) used by Borg [18] (square packing) is indicated by a dashed grey line and the theoretical MPD for mono-disperse circles in a rectangle is indicated by a solid grey line (hexagonal packing). The starting density for all simulations is indicated by the × symbol that is half visible on the x-axis.

Discussion
Our results show that the Hugoniot is directly affected by variations in grain arrangement (figures 9 and 11); as Z increases u p appears to decrease and U F increases. We use Spearman's rank correlation coefficient (represented by the symbol 'R S ' instead of ρ to avoid confusion with the symbol for density) to quantify these progressions. Spearman's rank assesses how well the relationship between the two variables is represented by a monotonic function [37,38]. R S can vary from −1, strong decreasing trend, to +1, strong increasing trend; a value of 0 implies no correlation. Both U F and u p have moderate to strong correlation with Z at all impact speeds. Correlation of u p with Z is strong for all V I ( R 0.9 S > | | , table 2) and the strength of correlation of U F against Z increases sharply between V I =450 ms −1 and V 550 ms Between these values it jumps from R S ≈0.5-0.6 (moderate) to R S ≈0.8-0.9 (strong), almost doubling (table 2). This suggests that coordination number Z is a useful predictor of the effect of mesostructure on shock response. However, shock velocity is intrinsically linked to the mechanisms of shock compaction as these dictate how the shock propagates, and so we examine whether this change in dependence of U F on Z, corresponds to a specific state of compaction. Figure 13 presents the Hugoniots in f−u p space, and in U F -u p space for comparison. Porosity was calculated in a similar manner to u p : f=1 − the mean of the volume fraction of cells in the shocked material. Borg and Vogler [18] proposed that the granular Hugoniot in the low u p range undergoes a transition (change in U s -u p slope) at a porosity of f≈21%, corresponding to the MPD for circles in a box, using square packing MPD=π/4≈0.79. MPD can also be defined using a hexagonal arrangement of circles, which gives a greater solid volume fraction of 0.91 The Hugoniot curves presented here display the two distinct slopes in the low u p region described in previous literature [18], but their transition does not occur as a sharp discontinuity at the square-packing MPD [18]. However, we note that the increase in correlation strength of U F against Z, from moderately correlated to strongly correlated, occurs at approximately the hexagonal packing MPD (f=9%, V I =450 ms −1 ). This observation suggests the possibility that when compaction exceeds the MPD the pre-shock Z becomes a primary predictor for the shock velocity and that the hexagonal MPD is a key state of compaction in the compaction of a particle bed. The effect of stiction on the shock response in numerical models has been discussed in previous literature, where removal of stiction in 2D simulations generally caused them to underestimate the shock response of an experiment for a given V I [11]. Here, we provide additional evidence for this. Our results show that removing stiction from the models results in an almost perfectly linear Hugoniot in U u F p space and the greater inter-particle motion generally produces post-shock densities greater than models with stiction. The inclusion of stiction produces a stiffer response (faster U F for a given u p ) and nonlinearities in the Hugoniot that are attributed to particle bed structure. Hence, for structure in the particle bed, characterised by Z, to influence the Hugoniot grain-grain contacts must provide some resistance to deformation.
In all the spaces considered, the least sensitive state variables to the mesostructure are longitudinal stress and particle velocity, displaying minimal absolute variation with Z. This insensitivity is well-known [36] and it can be inferred that the least sensitive plane to compare the shock response of particle beds with differing mesostructure will be the σ LL -u p plane ( figure 14). By minimising the effects of more sensitive variables it is possible to observe subtle changes in these insensitive variables. As V I increases the range of u p and LL s | | increase from ∼14 ms −1 and ∼0.03 GPa to ∼23 ms −1 and ∼0.1 GPa. However, the percentage range decreases from ∼13% and 8% to ∼4% and ∼2% respectively.
However, whilst u p and LL s | | have displayed the smallest absolute effect caused by changing Z, as previously discussed, U F and ρ have displayed the largest absolute effect. The effect of Z on the Hugoniots will be plainest in ρ-U F space as presented in figure 14. The trends are clearest in this space with the Hugoniots shifting to lower densities and higher shock velocities as Z increases, for constant V I .
The relative variation we observe in U s (ΔU s =2%-6%) and u p (Δu p =5%-15%), because of Z, is significant when compared to similar experiments. The percentage change in u p , Δu p , is greater than ΔU s because of the inherently smaller values involved; the absolute variation is smaller but we use relative variation here to compare directly with the experimental results. The original studies, that this study was based on, present experimental Hugoniot data with small uncertainties, but a large scatter. The uncertainty on U s was ≈0.5% [1], an order of magnitude less than the variation we observe (2%-6%). Additionally, for an approximately constant u p , the studies displayed a spread in U s of the order of ≈10% (figure 7) [27]. The scatter is of a comparable order of magnitude to the variation we observe due to Z but an order of magnitude larger than the reported uncertainty. This result suggests that the unquantified mesostructure could be a possible cause of this scatter. However, in studies of shock compaction of granular materials it is not common to quantify the mesostructure or grain arrangements of the pre-compacted samples, beyond measuring the bulk porosity. In experimental studies, it is practically impossible to use a common mesostructure for all experiments along a granular Hugoniot. Whilst it is possible to record the initial mesostructure in a particle bed before an experiment, it is typically difficult and costly: performing tomography on a sample before an experiment, for example. A simple method for dealing with the unknown variations in bed structure would be to combine results over multiple beds. In contrast, a record of pre-impact structure is significantly more practical to record in mesoscale simulation studies. However, in the case of simulations, this is typically not presented and algorithms for generating initial grain arrangements used in different studies are typically not the same and are sometimes not specified at all. The coordination number provides a simple metric to quantify mesostructure in a bed for comparison between studies. It is currently not feasible to calculate a value of Z in most experimental studies. However, a joint experimental-numerical approach could be employed in future studies to estimate the uncertainties caused by mesostructure variation.

Conclusions
We report a computational investigation into the effect of mesostructure on shock propagation through granular tungsten carbide. For the same bulk density, structure within the particle bed is defined by the number of contacts between particles or coordination number Z. The position of the Hugoniot in U F -u p space and LL s | |-ρ space is sensitive to the value of Z, and hence the amount of mesostructure in the particle bed, for the same initial porosity/bulk density. A bed with larger Z (more mesostructure) generally has a faster shock speed, lower post-shock bulk density, and lower particle velocity. For a given impact speed the particle and shock velocity in the moving frame (U F ) correlate strongly with Z. However the strength of correlation of U F with Z jumps from moderate to strong above an impact speed of 450 ms −1 . The transition corresponds approximately to the stress required to compact the bed to the equivalent hexagonal MPD for 2D circles in a box. Coordination number, therefore, appears to be a good predictor of shock response and its predictive efficacy increases significantly for stresses that result in near complete compaction of the particle bed. This may be because there is insufficient time for particle rearrangement at higher stresses. At lower stresses, involving low-tomoderate compaction, shock speed and bulk density appear to also depend on aspects of mesostructure beyond the simple measure of number of contact points. A more nuanced metric may be required to fully describe the dynamic compaction of a particle bed. Future work should consider the evolution of the coordination number with time within a bed; this is non-trivial because of the presence of deformation at many V I .
Our numerical results also show that the Hugoniot depends on stiction between grains. Stiction is an approximation of inter-particle shear resistance (friction and interlocking) that is used in Eulerian mesoscale simulations. The amount of stiction is varied by changing the strength of cells containing multiple materials, where two or more particles are in contact. Reducing this strength to zero, in all 'mixed' cells approximates removing inter-particle shear resistance entirely and produces a single linear U s -u p Hugoniot. The removal of stiction prevents particles locking together, making most existing contacts unstable structures and encourages particle motion into more stable arrangements. Different methods for treating strength in cells containing multiple materials will likely result in different Hugoniot curves, which should be accounted for when comparing results from mesoscale simulations.
This study has shown that varying the mesostructure within a mono-disperse particle bed has a distinct effect on its shock response. A simple method to account for this variation would be to combine results from multiple particle beds for the same model or experiment, but the coordination number, if tractable, provides a useful metric for estimating the effect. This work, however, focuses exclusively on the mono-disperse, circular, 2D scenario. Extending the investigation to 3D is necessary for developing a full understanding of the compaction mechanics. Studies have already shown that the transition is not trivial [11]. In addition, it is important to examine whether the effect of mesostructure on shock response observed here extends to poly-disperse, complex morphology particle beds. The coordination number may be an ideal metric for studying such cases but a greater understanding of its full effect, such as evolution in time, is required.