Simulations were done to gain insight whether the shock response of dry sand at low stresses would vary with porosity and whether the effects of friction between grains under confinement could be extracted from the planar plate impact experimental data. The sand sample was modeled as grains separated by voids representing porosity. The simulation procedure coupled grain deformations with frictional sliding at grain boundaries. The shock response of dry sand varied considerably with porosity. The sample compacted through pore closure followed by inelastic pore collapse mechanisms affecting the inhomogeneous response and shock rise time. The sample attained final compaction in the shock state long after attaining peak longitudinal velocity/stress. The calculated shock Hugoniot for a sample of high (40%) porosity was in agreement with experimental data. The Us-Up slopes for sand of 10% and 20% porosity were found to be negative. The calculated σH-ρH Hugoniot suggested that the two slopes would become positive at higher stresses in order to approach the solid Z-cut quartz Hugoniot at full compaction. High porosity sand may never exhibit negative slopes. It is concluded that the effects of friction between grains can be successfully extracted from a coupled experimental-computational approach. This requires measuring the velocity profile in the back buffer, elastic buffer material, and code capable of simulating frictional sliding between grains. The dispersion effect increased the slope of the velocity profile with propagation distance but did not result in a wave speed reduction or shock attenuation. This may be due to the small grain size and sample thickness as well as the absence of grain fragmentation in the present simulations.

In previous work,1 an ogive nose projectile was shown to become unstable and change direction when penetrating dry sand. This was due to inhomogeneous loading applied by the sand. The inherent heterogeneities of porous particulate sand lead to inhomogeneous loads as the projectile interacts with individual grains or clusters of grains during penetration. After initial impact, the impactor and secondary shock waves interact with sand that was compacted by the preceding waves. The medium resisting penetration ahead of the projectile showed spatial and random variation in the compacted density as observed in experiments.2–6 Compaction ahead of the projectile and its associated deformation field depends on the rate of penetration, the effects of friction,5,6 and the shape of the impactor nose.7,8 If the shock response of sand varies with porosity, then the resulting spatial variation in the sand response to shock loading and possible variation in frictional effects are possible sources of inhomogeneous loads on the projectile. Quantifying the shock response in penetration experiments is difficult and hence, attention is turned to the well-controlled planar plate impact experiments.

The shock responses of geologic materials including sand, ceramic, and concrete have been studied since the 1960s.9–19 The advent of mesoscale measurement techniques and multi-scale simulations20 has placed a renewed emphasis on the study of the shock response of dry and partially saturated sand21–32 and silica/ceramic powders.33–36 Free surface velocity profiles29 and manganin longitudinal stress gauge data18,30 elucidate their inhomogeneous response. Similarities between early9–19 and recent experiments21–32 include: (i) the initial porosity of sand is of the order of 40%, (ii) the measured shock velocity versus particle velocity (UsUp) or longitudinal stress versus particle velocity (σxUp) equation-of-states (EOS) are analogous to the EOS of conventional polycrystalline materials (e.g., metals), (iii) the slope of the linear UsUp relationship is positive, and (iv) the dispersion of the precursor wave30 is analogous to the attenuation of the elastic precursor in geologic materials.10 The minimum and maximum initial porosities achieved in recent studies are 33.8%29 and 95%,35 respectively. The in situ porosity of geologic materials, including sand, declines with depth and varies with location.37 A comparison of their EOS16 and porosity37 shows that the UsUp slopes of many geologic materials are negative at lower shock stresses in the 4%–20% porosity range and become positive at higher stresses. Similarly, many ceramic materials have a negative slope.13 Fused silica has a negative slope38,39 and a free volume (intrinsic porosity) in the same range.40 However, the shock response of pure sand with an initial porosity below 40% remains unknown due to the difficulties in preparing the sample through static compression.

One mechanism for the negative Us-Up slope is the transition to a quasi-hydrostatic state via dynamic yielding.16 Friction is known to be the predominant mechanism for yielding in geological media.41–43 However, the role of friction in the shock response of sand and whether the frictional relationship between grains under confinement can be extracted from measured data have not been addressed experimentally or computationally, to the best of our knowledge. Recent experimental studies are of limited use in this case because the inelastic deformation of aluminum, copper, and polymethyl-methacrylate (PMMA) sample encapsulating materials21–32 modifies the inhomogeneous response before it can be measured. The friction effect, if any, propagated by sand will be modified as well. For example, a PMMA backing plate measured24 a single shock with a sharp rise time at 0.983 km/s. In contrast, a copper backing plate21 measured an initial increase with a longer rise time followed by a dispersive increase at 0.969 km/s. The aluminum buffer disc27,28 presents similar difficulties with its inelastic deformation obscuring the dispersive rise time.28 Hence, to discern the friction effect, elastic encapsulating material is needed to ensure that the sand inhomogeneous response is propagated in its entirety to the measurement site.

Typical numerical modeling and simulation techniques have used homogeneous material models to simulate shock response. These include the Mie-Gruneisen equation-of-state,33 the Pα model,44 (and its variance26,27 implemented in CTH45 simulation code), the SP model35 (as a special case to the Pα model), the Pλ model,46 and pore collapse models.47–49 Stress is a function of porosity in these continuum models, and the two quantities simultaneously reach their final values in the peak shock state. When their parameters are determined from experimental data, the models predict the shock response of dry sand in reasonable agreement with data.26 They do not, however, predict dispersion in the velocity/stress profiles35 seen in the literature as oscillations during shock rise and in the peak shock state, increased shock rise time, and reduction in wave speed with propagation distance. This was demonstrated in earlier work12 that showed that even a detailed homogenized continuum model does not predict the dispersion and attenuation of the shock wave in concrete. The Pα and other void collapse models were developed to model shock compression in porous metals.50,51 These models assume homogeneous plastic deformation of the pore's cell walls as the mechanism for void collapse and porosity reduction. However, inelastic pore collapse in sand with an initial porosity as high as 40% can be envisioned only after enough compaction occurs to enable locking of the grains necessary for inelastic deformation. The limited mesoscale simulations of the past modeled individual metal powder particles as circular52 or rectangular.53 Recent mesoscale simulations27,36,54 also used circular particles to model sand or metal powder and predicted mesoscale response with varying success. Due to its ability to simulate large deformation of multi-material particulate systems, CTH45 has been the code of choice. However, due to the Eulerian formulation, CTH cannot model sliding with friction. It can model either no sliding (infinite friction) or free sliding (zero friction) at grain boundaries. This may be why the role of friction in the shock response of particulate media has not yet been computationally addressed.

The objectives of the present work were to simulate the shock response of porous dry sands in assumed plate impact experiments with a realistic model of mesoscale heterogeneities, determine whether the shock response at low stresses varied with initial porosity, and investigate whether the frictional relationship between grains under confinement could be extracted from experimental data. A brief description of the mesoscale simulation approach with response parameters, simulation results and discussions, and conclusions are given in Secs. IIIV.

Schematics of the simulated (a) plate impact experiment, (b) representative volume element (RVE), and (c) a zoomed mesoscale view of the sand sample with grains and pores are shown in Fig. 1. A 3 mm thick dry sand sample with grains and a given porosity was encapsulated between a 1 mm thick front cover plate (left buffer) and a 1 mm thick back plate (right buffer). The composite sample was backed by an 8 mm thick window and impacted by a 6 mm thick impactor. The impactor, buffer plates, and window were of homogeneous materials. The sand sample was generated by ISP-SAND55 software with an input mean grain size of 60 μm, ± 30% uniform size distribution, and an initial porosity up to 40%. Based on Voronoi tessellation,56 ISP-SAND generates grains of given mean size and size distribution mixed with a pre-calculated number of smaller grains (70% of the mean). The smaller grains, designated as pore grains, are removed to generate the initial porosity. The generated sand grains were 2D random polyhedrons with coincident edges. The zoomed mesoscale view of the sand with porosity is shown in Fig. 1(c). The coincident edges had no cohesion, and grains were allowed sliding with friction from the beginning of shock loading. The entire RVE including sand grains was discretized into constant strain triangles using CUBIT57 mesh generation software. Other parameters needed for simulations and required output parameters were specified using ISP-PREP software.58 

FIG. 1.

Schematic of the (a) assumed plate impact experiment, (b) zoomed view of the simulated RVE, and (c) mesoscale with grains in the sand sample. All units are in mm.

FIG. 1.

Schematic of the (a) assumed plate impact experiment, (b) zoomed view of the simulated RVE, and (c) mesoscale with grains in the sand sample. All units are in mm.

Close modal

The 2D updated Lagrangian explicit parallel finite element code, ISP2D,59 also referred to as ISP-TROTP,1,55,56 simulated the mesoscale response with lateral velocity in the Y-direction restrained at the two boundaries. All materials in the RVE, including sand grains, were assumed to be Z-cut quartz. This type of quartz is commonly used as window material27,60 in plate impact experiments and has well established isotropic shock Hugoniot properties.60 The abundantly available silica sand (SiO2) is composed mainly of quartz. Isotropic quartz data have been used to study the shock response of sand.3,29,32 The present assumption that the sand grains are isotropic Z-cut quartz grains is likely close to reality. The ISP2D code simulates the deformation of individual grains under shock loading using the specified elastic/inelastic materials model and grain-grain interaction using a contact algorithm with friction. The code does not yet have the capability to model grain fragmentation and thermal effects. The impactor, buffers, and window were elastic to enable the sand's inhomogeneous response to be propagated in its entirety to possible measurement sites. The elastic shock response was modeled by the non-linear Mie-Gruneisen equation-of-state for the volumetric response and mean stress dependent shear modulus for the deviatoric response.60 The same elastic models, coupled with an inelastic overstress model and a mean stress dependent yield strength, were used to simulate the shock response of the sand grains.61 The material models used in the simulations and their parameters are summarized in Table I. The impactor-buffer, grain-buffers, and grain-grain interactions were simulated using a contact algorithm1,55,56,59 with sliding and friction at their boundaries.

TABLE I.

Material models and their parameters used in simulations.

Z-cut quartz elastic model used for sand grains, impactor, buffer plates, and window 
Initial density: ρ0 = 2.6485 g/cm3 
Mean stress dependent shear modulus: G=46.92+1.873P+0.3459P2GPa46.92GPa 
P(=σm) is the negative of mean stress in GPa. 
Mie-Gruneisen equation of state: 
P=(43.19μ+156.2μ2+48.6μ3)(1γμ2)+γEV;μ=ρρ01;γρ=γ0ρ0 
μ is the volume compression, γ is the Mie-Gruneisen constant with γ0 = 0.675, E is the specific internal energy, and V is the specific volume. 
Inelastic overstress model for sand grains 
sij=sijT[ 1(3J2τy)ΔtTR ];J2=12sijsij; τy={ 2.353+0.667P0P3.067GPa4.4P>3.067GPaTR=3.5ns 
sij is the deviatoric stress, J2 is the second invariant of the deviatoric stress, τy is the mean stress dependent flow strength, Δt is the incremental time step, and TR is the relaxation time. 
Z-cut quartz elastic model used for sand grains, impactor, buffer plates, and window 
Initial density: ρ0 = 2.6485 g/cm3 
Mean stress dependent shear modulus: G=46.92+1.873P+0.3459P2GPa46.92GPa 
P(=σm) is the negative of mean stress in GPa. 
Mie-Gruneisen equation of state: 
P=(43.19μ+156.2μ2+48.6μ3)(1γμ2)+γEV;μ=ρρ01;γρ=γ0ρ0 
μ is the volume compression, γ is the Mie-Gruneisen constant with γ0 = 0.675, E is the specific internal energy, and V is the specific volume. 
Inelastic overstress model for sand grains 
sij=sijT[ 1(3J2τy)ΔtTR ];J2=12sijsij; τy={ 2.353+0.667P0P3.067GPa4.4P>3.067GPaTR=3.5ns 
sij is the deviatoric stress, J2 is the second invariant of the deviatoric stress, τy is the mean stress dependent flow strength, Δt is the incremental time step, and TR is the relaxation time. 

Simulations were carried out for 0%, 10%, 20%, 30%, and 40% input initial porosities of the sand sample. The 0% porosity sample of solid Z-cut quartz has a well-known60 shock response and is, therefore, not discussed. The ISP-SAND generated porosities that were within ±2.7% of the intended input porosities. Due to the absence of fragmentation and thermal effects, simulations were limited to a maximum impact velocity of 0.75 km/s, assuming that the two effects would be minimal and could be neglected. Sand samples with individual porosities were simulated at 0.25 km/s, 0.5 km/s, and 0.75 km/s impact velocities with a constant friction coefficient between grains of up to 0.8 with the initial mesoscale unchanged.

The shock response was quantified in terms of the line averaged temporal and area averaged spatial profiles of longitudinal velocity, stress, porosity, and density. The line averaged quantities were calculated from their values at 31 evenly spaced points on the 300 μm long tracer lines56 for the 60 μm mean grain size. As shown in Fig. 1(b), tracer lines L and R were at the center of the left and right buffers; I1 and I2 were at the sand interface with the left and right buffers, respectively; 1 to 5 were equally spaced in the sample. The area averaged longitudinal velocity [ =ΣVimi/Σmi ], longitudinal stress [ =(Σσxx(i)Ai)/A ], porosity [ =1(ΣAi)/A ], and density [ =(Σmi)/A ] were calculated by averaging over the rectangle of area A, where Ai is the area, mi is the mass, Vi is the longitudinal velocity, σxx(i) is the longitudinal stress at the centroid of triangle “I,” and summation was carried over triangles lying inside the averaging rectangle. The size of the rectangles and their positions used in the calculations are given with the results.

The shock Hugoniot was calculated in terms of the average peak particle velocity Up, shock velocity Us, density ρH, and longitudinal stress σH in the shock state. The shock velocity Us was calculated as Us=L0/(TI2TI1), where L0 is the initial sample thickness, and TI1 and TI2 are the shock arrival times at the left and right buffer interfaces, respectively. The arrival time was taken corresponding to half of the average peak velocity. The shock Hugoniot stress σH was calculated as σH=ρ0(1p0)UsUp. Density ρH in the shock state was calculated as ρH=ρ0(1p0)Us/(UsUp), where ρ0 is the solid Z-cut quartz density and p0 is the generated initial porosity. Unless otherwise noted, the simulation results discussed below and summarized in Table II were for samples of 60 μm mean grain size, ±30% uniform size distribution, elastic-inelastic grain deformation, and a 0.4 friction coefficient.

TABLE II.

Summary of simulation results and calculated shock Hugoniot parameters for dry sand with 60 μm mean grain size, ±30% uniform size distribution, elastic-inelastic grain response, and 0.4 friction coefficient between grain (unless otherwise mentioned).

Run No.Impact velocity (km/s)Up (km/s)Us (km/s)Peak ρ (g/cm3)ρH (g/cm3)Peak σxx (GPa)σH (GPa)
Input porosity 0%, solid Z-cut quartz 
0.25 0.125 6.500 2.700 2.700 2.151 2.152 
0.50 0.250 6.680 2.750 2.751 4.418 4.423 
0.75 0.389 6.690 2.820 2.812 6.552 6.892 
Input porosity 10%, generated porosity 9.8% 
0.25 0.141 5.690 2.462 2.449 1.880 1.917 
0.50 0.294 5.600 2.551 2.521 3.640 3.933 
0.75 0.460 5.080 2.676 2.627 5.190 5.583 
Input porosity 20%, generated porosity 19.5% 
0.25 0.172 4.290 2.252 2.222 1.368 1.574 
0.50 0.344 3.540 2.424 2.363 2.540 2.597 
0.75 0.540 3.230 2.568 2.561 3.700 3.720 
Input porosity 30%, generated porosity 29.8% 
10 0.25 0.215 1.660 2.180 2.137 0.667 0.664 
11 0.50 0.406 2.147 2.340 2.294 1.620 1.621 
12 0.75 0.599 2.441 2.470 2.465 2.700 2.720 
13a 0.50 0.416 1.974 2.392 2.357 1.376 1.527 
14b 0.50 0.406 2.214 2.320 2.278 1.679 1.672 
Input porosity 40%, generated porosity 40.0% actual porosity 
15 0.25 0.231 0.970 2.125 2.086 0.360 0.356 
16 0.50 0.443 1.580 2.280 2.208 1.112 1.112 
17 0.75 0.641 1.982 2.422 2.348 2.010 2.019 
Run No.Impact velocity (km/s)Up (km/s)Us (km/s)Peak ρ (g/cm3)ρH (g/cm3)Peak σxx (GPa)σH (GPa)
Input porosity 0%, solid Z-cut quartz 
0.25 0.125 6.500 2.700 2.700 2.151 2.152 
0.50 0.250 6.680 2.750 2.751 4.418 4.423 
0.75 0.389 6.690 2.820 2.812 6.552 6.892 
Input porosity 10%, generated porosity 9.8% 
0.25 0.141 5.690 2.462 2.449 1.880 1.917 
0.50 0.294 5.600 2.551 2.521 3.640 3.933 
0.75 0.460 5.080 2.676 2.627 5.190 5.583 
Input porosity 20%, generated porosity 19.5% 
0.25 0.172 4.290 2.252 2.222 1.368 1.574 
0.50 0.344 3.540 2.424 2.363 2.540 2.597 
0.75 0.540 3.230 2.568 2.561 3.700 3.720 
Input porosity 30%, generated porosity 29.8% 
10 0.25 0.215 1.660 2.180 2.137 0.667 0.664 
11 0.50 0.406 2.147 2.340 2.294 1.620 1.621 
12 0.75 0.599 2.441 2.470 2.465 2.700 2.720 
13a 0.50 0.416 1.974 2.392 2.357 1.376 1.527 
14b 0.50 0.406 2.214 2.320 2.278 1.679 1.672 
Input porosity 40%, generated porosity 40.0% actual porosity 
15 0.25 0.231 0.970 2.125 2.086 0.360 0.356 
16 0.50 0.443 1.580 2.280 2.208 1.112 1.112 
17 0.75 0.641 1.982 2.422 2.348 2.010 2.019 
a

0.1 Friction coefficient.

b

0.8 Friction coefficient.

Figure 2 shows the line averaged temporal profiles of (a) longitudinal velocity and (b) stress for 30% porous sand simulated at 0.5 km/s impact velocity. The inhomogeneous responses shown by oscillations during shock rise time and at peak shock state were comparable to those seen in recent data. The average peak shock state followed the one-dimensional (1D) shock propagation theory at a continuum scale. Stress was reduced and the velocity increased in the left buffer (L) by the release wave reflected from the interface (I1) with the low impedance sand sample. The peak velocity and stress in L and the sand (lines I1 to I3) were identical until the arrival of the recompression wave from interface (I2) with a higher impedance buffer (R). Ultimately, the average peak velocity and stress at I2 and R were identical.

FIG. 2.

Line averaged in situ temporal profiles of (a) longitudinal velocity and (b) longitudinal stress at 0.5 km/s impact velocity.

FIG. 2.

Line averaged in situ temporal profiles of (a) longitudinal velocity and (b) longitudinal stress at 0.5 km/s impact velocity.

Close modal

The two vertical lines in Fig. 2(a) show the time when the peak shock state is reached at lines 1 and 3 (0.677 μs and 1.402 μs, respectively) in terms of longitudinal velocity, well before the arrival of the recompression wave from the right buffer. Features resembling precursors at the foot of the velocity profiles are shown in Fig. 2(a). These precursor-like features were present in all simulation results and varied randomly when the sand mesoscale was changed through ISP-SAND generation (grain size and porosity were unchanged). The feature disappeared at some tracer lines and its magnitude varied randomly along the shock propagation direction. Numerical artifacts associated with the use of artificial viscosity59 in simulations also contributed to the slow increase in velocity and stress, yielding a precursor like feature. The small sample thickness, small grain size, randomness, and numerical artifacts prevented us from identifying the feature as a precursor or determining its dispersion and attenuation. The shock rise time in the sand was higher than that of conventional polycrystalline materials (metals and ceramics) and increased with propagation distance due to dispersion. Figure 3 shows the dispersion effect through a more apparent increase in slope of the area averaged spatial velocity profile with propagation distance. The profiles were calculated by averaging over a 50 μm thick rectangle and sweeping the sample along the longitudinal x-direction at discrete times. However, no discernable reduction in wave speed was found. Hence, calculation of the average peak particle velocity from the in situ line averaged velocity profiles (lines 1–3) and shock velocity from the time it took for the wave to travel from interface I1 to I2 was justified. The shock Hugoniot parameters calculated from the simulation were: Up = 0.406 km/s, Us = 2.147 km/s, σH = 1.62 GPa, and ρH = 2.34 g/cm3. The average peak shock stress (lines 1–3) remained constant, showing a lack of attenuation. The constant wave speed and lack of attenuation were further confirmed (not shown) by simulating a 6 mm thick sample with 120 μm average grain size (other parameters were unchanged). The average peak velocity and stress in the right buffer were the same as that of the 3 mm thick sample.

FIG. 3.

Area averaged in situ spatial profiles of longitudinal velocity at discrete times simulated at 0.5 km/s impact velocity.

FIG. 3.

Area averaged in situ spatial profiles of longitudinal velocity at discrete times simulated at 0.5 km/s impact velocity.

Close modal

Figure 4 shows the deformed mesoscale with longitudinal stress contours at 1.5 μs. Two mechanisms of compaction were apparent in the figure and were found in all simulations independent of porosity. After entering the sample, the shock wave loaded the grains in compression adjacent to interface I1. The compressed grains unloaded partly by reflection from their free surfaces adjacent to the pores and gained momentum impinging on the next grain. This led to pore closure as the first mechanism signified by grain translations in both directions without significant energy loss and predominantly elastic grain deformations. The pore closure mechanism continued until the grains were locked by friction forming a force chain. This resulted in inelastic pore collapse as the second mechanism signified by inelastic grain deformations, relative sliding, and rotations. The elastic pore closure mechanism was predominant at higher porosity and lower impact velocity. In contrast, the inelastic pore collapse mechanism was predominant at lower porosity and higher impact velocity. The spatially distributed inelastic compaction slowed the rate of velocity increase before reaching the peak state causing a significant increase in the rise time as evident in Fig. 2(a). The shock rise time decreased with the increase in impact velocity (not shown). In other words, dispersion was more profound at higher porosities and lower impact velocities.

FIG. 4.

Deformed mesoscale of sand sample showing force chains and longitudinal stress contours at 1.5 μs simulated at 0.5 km/s impact velocity.

FIG. 4.

Deformed mesoscale of sand sample showing force chains and longitudinal stress contours at 1.5 μs simulated at 0.5 km/s impact velocity.

Close modal

Figure 5 shows the in situ area averaged porosity and density calculated by averaging over 300 μm thick rectangles centered at tracer lines 1, 3, and 5. The two vertical lines (Fig. 2(a)) show the times of attaining peak shock state at tracer lines 1 and 3. The rate of compaction was higher during the initial pore closure and decreased during inelastic pore collapse as seen for lines 1 and 3. A comparison of Fig. 5 with Fig. 2(a) shows that shock rise time was shorter during pore closure and compaction continued long after the peak longitudinal velocity/stress state was reached. The first two areas at tracer lines 1 and 3 had similar final porosity of 13.5% and had average densities the same as the calculated ρH in the shock state. The recompression wave arrived at tracer line 5 before the in situ shock state was reached (Fig. 2). As a result, the initial compaction under the combined effect of the two waves was longer and final porosity was lower. However, the rate of compaction declined before reaching the final porosity signifying increased inelastic deformation. The elastic only grains (by suppressing inelasticity) compacted by pore closure with shorter rise time, increased inhomogeneous response, and reduced compaction (not shown). Thus, inelastic compaction was the primary reason for slow rise time and rounding of velocity profiles.

FIG. 5.

Area averaged in situ temporal profiles of the porosity and density at tracer lines 1, 3, and 5 (left to right) at 0.5 km/s impact velocity.

FIG. 5.

Area averaged in situ temporal profiles of the porosity and density at tracer lines 1, 3, and 5 (left to right) at 0.5 km/s impact velocity.

Close modal

Figure 6 shows the line averaged longitudinal velocity profiles at L and R simulated at 0.5 km/s impact velocity with 10%, 20%, and 40% initial porosities. The results for the 30% porosity discussed above are also shown for comparison. The increase in porosity reduced the sample density and shock impedance, resulting in wave speed reduction, higher peak particle velocity, lower in situ stress (not shown), lower peak particle velocity in R, and increased inhomogeneous response. The 10% porous sand compacted primarily through an inelastic mechanism with a sharp rise time and fewer oscillations. The contribution of pore closure increased in 20% porous sand causing a gradual increase in velocity. The 30% porous sand compacted through pore closure first, with a sharp rise time followed by the inelastic process discussed above. Pore closure was the primary mechanism in the 40% porous sand and resulted in a predominantly elastic grain response, a shorter rise time, and greater inhomogeneous response and dispersion.

FIG. 6.

Line averaged in situ temporal profiles of longitudinal velocity with different initial porosity simulated at 0.5 km/s impact velocity.

FIG. 6.

Line averaged in situ temporal profiles of longitudinal velocity with different initial porosity simulated at 0.5 km/s impact velocity.

Close modal

The sand samples with different porosities were simulated at 0.25 km/s, 0.5 km/s, and 0.75 km/s impact velocities keeping the initial mesoscale geometry unchanged. The calculated shock Hugoniot parameters are summarized in Table II. Figure 7 shows the calculated (a) UsUp and (b) σHρH shock Hugoniots for the five initial porosities. The UsUp relations at different porosities are

US=6.443+0.707UPfor   p0=0%,Up0.45km/s,US=6.0321.929UPfor   p0=10%,Up0.48km/s,US=4.6902.850UPfor   p0=20%,Up0.56km/s,US=1.256+2.033UPfor   p0=30%,Up0.62km/s,US=0.427+2.472UPfor   p0=40%,Up0.65km/s.
(1)
FIG. 7.

Shock Hugoniot parameters (a) Us-Up, and (b) σH-ρH calculated from the mesoscale simulation results.

FIG. 7.

Shock Hugoniot parameters (a) Us-Up, and (b) σH-ρH calculated from the mesoscale simulation results.

Close modal

Fig. 7(a) also shows the literature data US=0.243+2.348UP for 41% initial porosity26 and US=0.511+1.72UP for 45% initial porosity.21 The calculated UsUp for 40% porosity agreed well with the literature data. The shock response of dry sand at low stresses significantly depends on initial porosity. The simulations predicted a negative slope for the UsUp relation in 10% and 20% porous sand for the first time. The grain deformations for the two porosities were predominantly elastic at low velocity and transitioned to predominantly inelastic at higher impact velocities. This suggests a possible mechanism for the negative slope—a transition from a faster elastic wave to the slower inelastic wave reducing Us at higher Up. A clear elastic precursor wave followed by an inelastic wave was not predicted due to the inability of porous media to support and propagate single shock waves. The solid Z-cut quartz (0% porosity) has an elastic response at low stress and a positive UsUp slope. The UsUp slope is positive for 30% and 40% porosities, and the sand response was largely elastic at these two porosities. The calculated σHUp relationship (not shown) indicated that the Hugoniot stress σH was on the release paths for the solid Z-cut quartz within ±7% and agreed with experimental data.32 The right-most line in Fig. 7(b) is the σHρH plot for the solid Z-cut quartz. The gradual concave curves for 10% and 20% porosities represent the transition from elastic to inelastic deformations and a negative UsUp slope. If it is assumed that the sand sample has to compress with increasing shock stress and approach the solid line for Z-cut quartz at full compaction, then the slope of the curves for the 10% and 20% porosities will become convex to approach the solid line. Hence, the UsUp slope for the two porosities will change from negative to positive before approaching the fully compressed state. The shock stress needed to fully compress the 40% porous sand will be high enough for the inelastic waves to be overdriven keeping the UsUp slope positive. Hence, the high porosity sand may never show a negative UsUp slope.

The line averaged (a) longitudinal velocity and (b) stress profiles in the left and right buffers simulated at 0.5 km/s impact velocity with 0.1 and 0.8 friction coefficients are shown in Fig. 8. For comparison, previous results for a friction coefficient of 0.4 are also shown. The initial mesoscale geometry with 30% porosity was the same for the three simulations. The average peak longitudinal velocity and stress were the same within tolerance for the three friction values. The wave speed and shock rise time were lower at the low friction coefficient of 0.1 due to the increased inhomogeneous response, dispersion, reduced inelastic deformation, and the increased pore closure mechanism of compaction. The compaction was more uniform with reduced porosity (≈12%) in the peak shock state. The increase in friction accelerated grain locking and increased the inelastic response. The early setup of a force chain provided medium continuity for early wave arrival at the right buffer. Less uniform compaction due to the smaller lateral motion resulted in higher porosity (≈15.5%) in the peak shock state. The longer shock rise time with greater friction increased the slope of the loading portion of the velocity and stress profiles. The distinct increase in the profile slopes suggests that the friction value between grains under confinement can be extracted from a combined experimental-simulation approach as long as the velocity profile is measured in the right buffer, the buffer is elastic, and the simulation code can simulate sliding with friction between grains. A series of simulations with different friction values will yield the average friction between grains that predicts the velocity profile in the right buffer in agreement with experimental data.

FIG. 8.

Line averaged temporal profiles of (a) longitudinal velocity and (b) longitudinal stress with different friction values between grains simulated at 0.5 km/s impact velocity.

FIG. 8.

Line averaged temporal profiles of (a) longitudinal velocity and (b) longitudinal stress with different friction values between grains simulated at 0.5 km/s impact velocity.

Close modal

Two-dimensional mesoscale simulations were carried out to gain insight into the shock response of dry sand at low stresses with varying initial porosity and friction values. The simulations of a planar plate impact experiment used a realistic model of the heterogeneities of sand sample. The elastic and elastic-inelastic deformations of sand grains coupled with friction and sliding at the grain boundaries, were modeled using ISP2D. The impactor, sample encapsulating buffers, and window were input as elastic materials to enable the sand's inhomogeneous response to be propagated in its entirety. As a result, the simulated response was more inhomogeneous than previously reported in the literature.

Oscillations in longitudinal velocity/stress during the rise and at peak shock state and elongated shock rise time represented a typical inhomogeneous response. The peak average in situ velocity and stress followed the general one-dimensional continuum shock propagation theory. The inhomogeneous response was strongly dependent on porosity, impact velocity, and friction between grains. The effect of grain size was not presented because its variation was coupled with variations in pore size. The ISP-SAND program generates a pre-calculated number of smaller pore grains, which are removed to generate pores. There is a trade-off in assuming pore grain size. Too small a fraction of mean grain size increases the number of pore grains for the input porosity. It becomes difficult to seed the RVE with more pore grain sites, causing the generated porosity to be well below its input value. Too large a fraction gives the desired porosity, but leads to unrealistically large pores. A realistic pore distribution resulted from a pore grain size of 70% of the mean grain size. As a result, the grain size effects were always confounded with pore size effects.

Two distinct compaction mechanisms were found in the simulations. Sand compacted by pore closure mechanisms first with grain translations, impingement, and primarily elastic grain deformations. This was followed by an inelastic pore collapse mechanism with inelastic deformations and rotations of grains after the frictional forces locked the grains together to form force chains. The two mechanisms were strongly dependent on porosity, impact velocity, and friction value. The pore closure mechanism was predominant at higher porosities and lower impact velocities, while inelastic pore collapse was predominant at lower porosities and higher impact velocities. Their coupled effect significantly increased the shock rise time. The rate of compaction by pore closure was greater than that of the inelastic collapse mechanism. More importantly, the compaction process continued long after the average peak longitudinal velocity and stress were reached. This means that the evidence of average peak velocity (or stress) does not necessarily mean that the sand had attained a true shock state, making the sample thickness an important parameter. Simulations with sufficiently thick sample will allow determination of time sand characteristically takes to reach the final compacted density after attaining the peak shock state. This characteristic time may depend on initial grain size and fragmentation during shock loading.

The shock response of dry sand at low stresses varied considerably with initial porosity. The calculated shock Hugoniot had a negative UsUp slope at low porosities and a positive slope at higher porosities. For the first time, the simulations predicted negative UsUp slopes for 10% and 20% initial porosities—analogous to other geologic media with the same porosity range. The calculated σHρH relationships showed that the UsUp slope at 10% and 20% porosities would need to become positive at higher stresses for the fully compacted state to approach the solid Z-cut quartz Hugoniot. However, the 40% porous sand would never have a negative slope due to the overdriven inelastic wave at very high stresses needed for full compaction. The simulations were not carried out at higher impact velocities to confirm the inferred change in slope due to unrealistic and excessively large inelastic deformation of grains, which created instabilities in the contact algorithm and simulation breakdown. Grain fragmentation is the real phenomena and must be modeled together with the thermal effect in simulations at higher impact velocities.

We conclude that friction between grains under confinement can be extracted from a combined experimental-computational approach if the velocity is measured in the right buffer. Accelerated grain locking and increased inelastic deformation at high friction values appreciably increased the shock rise time and slope of the velocity profile. This suggests the need for elastic sample encapsulating material; otherwise, the inelastic deformation obscures the above effects. The simulation code to extract the friction values should be capable of modeling sliding with friction between grains. Moreover, any future plate impact experiments to find the friction values at different impact velocities or repeat experiments should use single-source sand to prevent variation in grain surface characteristics. Samples should be prepared with identical grain size distributions, compaction processes, and porosities to prevent excessive variation in frictional effects from one experiment to the next.

The shock wave dispersion effect in sand was limited to an increase in the slope of the velocity/stress profiles with propagation distance. The simulations did not show a reduction in wave speed or shock attenuation as found in other geologic materials and concrete. The reasons may be small grain size, small sample thickness, and lack of grain fragmentation effects in the simulation procedure. Grain fragmentation has yet to be modeled, to the best of our knowledge. The actual size of sand grains varies from 80 μm to more than 1 mm. Such large grain sizes would require a sufficiently thick sample to ensure that the response is measured or simulated for shock propagation through a sufficient number of grains in the longitudinal direction to provide an averaged response. Furthermore, the sample thickness has to provide a sufficient propagation length to study dispersion and attenuation before the arrival of the recompression wave from the right buffer. Such large simulations were not possible due to computational resource requirements. Finally, variation in the shock response of dry sand due to porosity and possible variations in friction due to the spatial distribution of compaction may be key contributing factors to the inhomogeneous projectile loads. Plate impact experimental data below 33.8% initial porosity and evidence of a negative slope for dry sand at lower porosities do not exist in the literature; this may be due to difficulties in sample compaction.

In the end, limitations of the present work need to be mentioned to drive the importance of future work. First, the present work is based on 2D plane strain simulations, while the shock response of heterogeneous particulate sand is essentially three-dimensional (3D) and requires 3D simulations. Three-dimensional simulations require greater computational resources by an order of magnitude and may yield new shock response features not captured by 2D simulation. However, 1D shock propagation theory is used to extract the shock Hugoniot of porous sand from plate impact experimental data. The present 2D simulations predicted the shock Hugoniot in reasonable agreement with experimental data. This indicates that the average response at continuum scale derived from 2D and 3D mesoscale simulations should agree with one another and with the experimental data. Second, the sand grains generated and used in these simulations were polyhedrons with coincident edges, which may not be realistic. The polyhedron closely resembles the shape of the sand grain. However, the coincident edges were the result of a Voronoi tessellation algorithm. The coincident edges did not have cohesion and grains were free to move, rotate, and impinge on other grains. A new algorithm is under development that will eliminate the need for pore size, thus decoupling pore size from grain size. The new algorithm will generate grains in random contact with one another for the given grain size, size distribution, and porosity. Another approach is to use real mesoscale measured by X-ray tomography62,63 in the simulations. While novel, the approach requires an actual sand sample, limiting variation in individual mesoscale parameters (grain size, size distribution, pore size, porosity, etc.). While any computational approach should ultimately simulate the sample with real mesoscale, adopting real mesoscale as the sole approach will limit the computational research. For example, experiment for sand with sample initial porosity less than 33.8% has not been possible due to the difficulties in static compaction process of preparing sample. But, computational process generated sample with initial porosities of 10% and 20% and negative UsUp slope was predicted for the two porosities for the first time. Hence, simulations with computationally generated sample as well as real sample from X-ray tomography may be essential to aid experiments as well as simulations. Finally, and more importantly, to the best of our knowledge, grain fragmentation and thermal effects during shock loading of sand have not yet been modeled. The use of small grain size and low impact velocities in simulations offset the associated disadvantages to some extent. However, there is a need for a robust computational approach to model grain fragmentation and thermal effects during shock loading of particulate media. Only then will realistic simulations be possible at higher impact velocities leading to full compaction of sand.

This work was supported by a DOE/NNSA grant. Professor Y. M. Gupta is acknowledged for suggesting this work, his support and many insights, and numerous useful discussions during the course of this work.

1.
S. K.
Dwivedi
,
R.
Teeter
,
C.
Felice
, and
Y. M.
Gupta
,
J. Appl. Phys.
104
,
083502
(
2008
).
2.
K.
Watanabe
,
K.
Tanaka
,
K.
Iwane
,
S.
Fukuma
,
K.
Takayama
, and
H.
Kobayashi
,
Technical Report No. AOARD-094011
, Chubu University, Japan,
2010
.
3.
J. P.
Borg
,
M. P.
Morrissey
,
C. A.
Perich
,
T. J.
Vogler
, and
L. C.
Chhabildas
,
Int. J. Impact Eng.
51
,
23
(
2013
).
4.
J. W.
Addiss
,
A.
Collins
,
F.
Bobaru
,
K.
Promratana
, and
W. G.
Proud
, in
9th International Conference on the Mechanical and Physical Behaviour of Materials under Dynamic Loading
,
Brussels
,
Belgium
,
2009
, Vol.
1
.
5.
J. W.
Addiss
,
A.
Collins
, and
W. G.
Proud
, in
Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Nashville, TN
, edited by
M. L.
Elert
,
W. T.
Buttler
,
M. D.
Furnish
,
W. W.
Anderson
, and
W. G.
Proud
(
2009
), p.
1313
.
6.
J. W.
Addiss
,
A.
Collins
, and
W. G.
Proud
, in
Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Nashville, TN
, edited by
M. L.
Elert
,
W. T.
Buttler
,
M. D.
Furnish
,
W. W.
Anderson
, and
W. G.
Proud
(
2009
), p.
1357
.
7.
A.
Collins
,
J.
Addiss
, and
W.
Proud
, in
Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Nashville, TN
, edited by
M. L.
Elert
,
W. T.
Buttler
,
M. D.
Furnish
,
W. W.
Anderson
, and
W. G.
Proud
(
2009
), p.
1443
.
8.
A. L.
Collins
,
J. W.
Addiss
,
S. M.
Walley
,
K.
Promratana
,
F.
Bobaru
,
W. G.
Proud
, and
D. M.
Williamson
,
Int. J. Impact Eng.
38
,
951
(
2011
).
9.
T. J.
Ahrens
and
V. G.
Gregson
,
J. Geophys. Res.
69
,
4839
, doi: (
1964
).
10.
T. J.
Ahrens
and
G. E.
Duvall
,
J. Geophys. Res.
71
,
4349
, doi: (
1966
).
11.
“Compendium of shock wave data,” UCRL-50108, Physics TID-4500, edited by
M.
van Theil
,
A. S.
Kusubov
, and
A. C.
Mitchell
, Vol. 1, Sec. A-2, Lawrence Radiation Laboratory,
University of California
,
Livermore, USA
,
1977
.
12.
Y. M.
Gupta
and
L.
Seaman
,
Nucl. Eng. Des.
45
,
507
(
1978
).
13.
LASL Shock Hugoniot Data
, edited by
S. P.
Marsh
(
University of California Press
,
1980
).
14.
B. M.
Lempriere
, “
Critical mechanics problems in system design
,” in
Proceedings of the Army Symposium on Solid Mechanics
(
Bass River
,
MA, USA
,
1982
).
15.
J. K.
Gran
,
Y. M.
Gupta
, and
A. L.
Florence
,
Mech. Mater.
6
,
113
(
1987
).
16.
T. J.
Ahrens
and
M. L.
Johnson
,
Rock Physics and Phase Relation: A Handbook of Physical Constants
(
The American Geophysical Union
,
1995
), Vol.
35
.
17.
G.
Yuan
,
R.
Feng
, and
Y. M.
Gupta
,
J. Appl. Phys.
89
,
5372
(
2001
).
18.
K.
Tsembelis
,
D. J.
Chapman
,
C. H.
Braithwaite
, and
W. G.
Proud
, in
Proceedings of the 2006 SEM Annual Conference and Exposition
,
2006
.
19.
G. R.
Willmot
and
W. G.
Proud
,
Int. J. Rock Mech. Min. Sci.
44
,
228
(
2007
).
20.
M. R.
Baer
and
W. M.
Trott
, in
Proceedings of the 12th International Symposium on Detonation
(
Office of Naval Research
,
2005
), Vol.
333-05-2
.
21.
K.
Tsembelis
,
W. G.
Proud
,
B. A. M.
Vaughan
, and
J. E.
Field
, in
Proceedings of the 13th DYMAT Technical Meeting, Saint-Jean-de-Luz, France
(
2001
), p.
193
.
22.
A. D.
Resbyansky
and
N. K.
Bourne
, in
Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter
, edited by
M. D.
Furnish
,
Y. M.
Gupta
, and
J. W.
Forbes
(
American Institute of Physics
,
Portland, Oregon
,
2003
), p.
1474
.
23.
A. D.
Rsenyansky
and
N. K.
Bourne
,
J. Appl. Phys.
95
,
1760
(
2004
).
24.
D. J.
Chapman
,
K.
Tsembelis
, and
W. G.
Proud
, in
Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Baltimore, MD, USA
, edited by
M. D.
Furnish
,
M.
Elert
,
T. P.
Russell
, and
C. T.
White
(
2005
), p.
1445
.
25.
D. J.
Chapman
,
K.
Tsembelis
, and
W. G.
Proud
, in
Proceedings of the 2006 SEM Annual Conference and Exposition
,
2006
.
26.
J. L.
Brown
,
T. J.
Vogler
,
D. E.
Grady
,
W. D.
Reinhart
,
L. C.
Chhabildas
, and
T. F.
Thornhill
, in
Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Waikoloa, Hawaii, USA
, edited by
M.
Elert
,
M. D.
Furnish
,
R.
Chau
,
N.
Holmes
, and
J.
Nguyen
(
2007
), p.
1363
.
27.
J. L.
Brown
,
T. J.
Vogler
,
L. C.
Chhabildas
,
W. D.
Reinhart
, and
T. F.
Thornhill
, Sandia National Laboratory Report No. SAND2007-3524,
2007
.
28.
T. J.
Vogler
,
M. Y.
Lee
, and
D. E.
Grady
,
Int. J. Solids Struct.
44
,
636
(
2007
).
29.
M.
Arlery
,
M.
Gardou
,
J. M.
Fleureau
, and
C.
Mariotti
,
Int. J. Impact Eng.
37
,
1
(
2010
).
30.
W. D.
Neal
,
D. J.
Chapman
, and
W. G.
Proud
,
Eur. Phys. J.: Appl. Phys.
57
,
31001
(
2012
).
31.
N.
Zdurr
,
M.
Sauer
,
N.
Guldemeister
,
K.
Wunnemann
, and
S.
Heiermaier
, in
Proceedings of the 12th Hypervelocity Impact Symposium
,
Baltimore
,
MD, USA
, edited by
B.
Sorenson
,
2012
.
32.
C. H.
Braithwaite
,
J. I.
Perry
,
N. E.
Taylor
, and
A. P.
Jardine
,
Appl. Phys. Lett.
103
,
154103
(
2013
).
33.
J. P.
Borg
,
L.
Schwalbe
,
J.
Cogar
,
D. J.
Chapman
,
K.
Tsembelis
,
A.
Ward
, and
A.
Lloyd
, in
Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Baltimore, MD, USA
, edited by
M. D.
Furnish
,
M.
Elert
,
T. P.
Russell
, and
C. T.
White
(
2005
), p.
37
.
34.
J. P.
Borg
,
D. J.
Chapman
,
K.
Tsembelis
,
W. G.
Proud
, and
J. R.
Cogar
,
J. Appl. Phys.
98
,
073509
(
2005
).
35.
J. P.
Borg
,
J. R.
Cogar
,
A.
Lloyd
,
A.
Ward
,
D. J.
Chapman
,
K.
Tsembelis
, and
W. G.
Proud
,
Int. J. Impact Eng.
33
,
109
(
2006
).
36.
J. P.
Borg
and
T. J.
Vogler
,
Int. J. Solids Struct.
45
,
1676
(
2008
).
37.
G. E.
Manger
, Geological Survey Bulletin, 1144-E, U.S. Atomic Energy Commission,
1963
.
38.
R.
Feng
and
Y. M.
Gupta
,
Materials Models for Sapphire, Alpha-Quartz, Lithium Fluoride, and Fused Silica for Use in Shock Wave Experiments and Wave Code Calculations
(
Institute for Shock Physics Internal Report, Washington State University
,
Pullman, WA
,
1996
).
39.
C. S.
Alexander
,
L. C.
Chhabildas
,
W. D.
Reinhart
, and
D. W.
Templeton
,
Int. J. Impact Eng.
35
,
1376
(
2008
).
40.
R. H.
Doremus
,
J. Am. Ceram. Soc.
49
,
461
(
1966
).
41.
T.
Ramamurthy
,
Int. J. Rock Mech. Min. Sci.
38
,
683
(
2001
).
42.
C. A.
Bareither
,
T. B.
Edil
,
C. H.
Benson
, and
D. M.
Mickelson
,
J. Geotech. Geoenviron. Eng.
134
,
1476
(
2008
).
43.
R.
Alias
,
A.
Kasa
, and
M. R.
Taha
,
Int. J. Civil, Architectural, Struct. Construction Eng.
8
,
1084
(
2014
).
44.
W.
Hermann
,
J. Appl. Phys.
40
,
2490
(
1969
).
45.
R. L.
Bell
,
M. R.
Baer
,
R. M.
Brannen
,
M. G.
Elrich
,
E. S.
Hertel
,
S. A.
Silling
, and
P. A.
Taylor
, “
CTH user's manual and input instructions
,” Sandia National Laboratory Report,
2003
.
46.
D. E.
Grady
,
N. A.
Winfree
,
G. I.
Kerley
,
L. T.
Wilson
, and
L. D.
Kuhns
,
J. Phys. IV
10
,
15
(
2000
).
47.
M. M.
Caroll
and
A. C.
Holt
,
J. Appl. Phys.
43
,
1626
(
1972
).
48.
V. F.
Nesterenko
,
Dynamics of Heterogeneous Materials
(
Springer-Verlag
,
New York
,
2001
).
49.
W.
Tong
and
G.
Ravichandran
, in
High Pressure Shock Compression of Solids: Response of Highly Porous Solids to Shock Loading
, edited by
L.
Davison
,
Y.
Horie
, and
M.
Shahinpur
(
Springer-Verlag
,
New York
,
1997
), p.
177
.
50.
B. M.
Butcher
and
C. H.
Karnes
,
J. Appl. Phys.
40
,
2967
(
1969
).
51.
J. R.
Asay
,
J. Appl. Phys.
46
,
197
(
1975
).
52.
J. E.
Flinn
,
R. L.
Williamson
,
R. A.
Berry
,
R. N.
Wright
,
Y. M.
Gupta
, and
M.
Williams
,
J. Appl. Phys.
64
,
1446
(
1988
).
54.
T. J.
Vogler
,
C. S.
Alexander
,
J. L.
Wise
, and
S. T.
Montgomery
,
J. Appl. Phys.
107
,
043520
(
2010
).
55.
R.
Teeter
, MS Thesis, “
Two dimensional mesoscale simulations of projectile instability during penetration of dry sand
,”
Washington State University
,
2007
.
56.
S. K.
Dwivedi
,
J. R.
Asay
, and
Y. M.
Gupta
,
J. Appl. Phys.
100
,
083509
(
2006
).
57.
"CUBIT Mesh Generation Environment Volume 1: Users Manual," Sandia National Laboratory Report No. SAND94-1100.
58.
S. K.
Dwivedi
,
ISP-PREP: Preprocessor for the 2D Finite Element Simulation of the Plate Impact Material Property Experiment
(
Institute for Shock Physics, Washington State University
,
Pullma, WA
,
2007
).
59.
S. K.
Dwivedi
and
Y. M.
Gupta
, “
ISP2D: Two-dimensional Lagrangian finite element code to simulate shock wave propagation in homogeneous and heterogeneous media
,”
Institute for Shock Physics Report No. ISP-TR-09-1
,
Washington State University
,
Pullman, WA
,
2009
.
60.
J. M.
Winey
,
R.
Feng
, and
Y. M.
Gupta
, “
Isotropic material models for the elastic response of sapphire and quartz single crystals under shock wave loading
,”
Institute for Shock Physics Internal Report No. 01-XX
,
Washington State University
,
Pullman, WA
,
2001
.
61.
C. H. M.
Simha
and
Y. M.
Gupta
,
J. Appl. Phys.
96
,
1880
(
2004
).
62.
C.
Halverson
,
D. J.
White
, and
J.
Gray
, in
The 2005 Mid-Continent Transportation Research Symposium
,
Ames
,
Iowa
,
2005
.
63.
B.
Zhou
,
J.
Wang
, and
B.
Zhao
,
Eng. Geol.
184
,
126
(
2015
).