Development and validation of model for sand

There is a growing requirement within QinetiQ to develop models for assessments when there is very little experimental data. A theoretical approach to developing equations of state for geological materials has been developed using Quantitative Structure Property Modelling based on the Porter-Gould model approach. This has been applied to well-controlled sand with different moisture contents and particle shapes. The Porter-Gould model describes an elastic response and gives good agreement at high impact pressures with experiment indicating that the response under these conditions is dominated by the molecular response. However at lower pressures the compaction behaviour is dominated by a micro-mechanical response which drives the need for additional theoretical tools and experiments to separate the volumetric and shear compaction behaviour. The constitutive response is fitted to existing triaxial cell data and Quasi-Static (QS) compaction data. This data is then used to construct a model in the hydrocode. The model shows great promise in predicting plate impact, Hopkinson bar, fragment penetration and residual velocity of fragments through a finite thickness of sand.


Introduction
QinetiQ has a growing requirement for developing assessments of granular materials when there is little experimental data to fit a model.This has necessitated the development of theoretical tools which have universal application to a wide range of materials.
This paper describes the use of Quantitative Structure-Property Modelling (QSPM) in deriving physically based equations of state directly, from knowledge of the constituents of the material.For geological materials this is based on an implicit assumption that many geological materials are a derivative of different crystal forms of silica (i.e.crystobalite, coesite, α-quartz, stishovite).The constitutive model is then fitted to existing Quasi-Static (QS) triaxial cell data.The ensuing model is then tested against a range of experiments covering different regimes of behaviour.The main issue with this approach is the compaction behaviour which is very variable and difficult to quantify.

Quantitative structure property modelling of silica
Silica is a molecule which exhibits 4 main crystal forms.These are in ascending order of density: crystobalite, quartz, coesite and stishovite and their molecular structure is illustrated in Fig. 1.CASTEP [1] reproduced the minimum energy structures of the different polymorphs of silica with the specified space group of the crystal forms and predicted their correct density values.However, CASTEP was unable to simulate changes in the crystal a Corresponding author: pdchurch@qinetiq.comstructure of silica with different boundary conditions such as pressure or volume imposed.The practical consequence of this problem was that the predicted bulk elastic modulus of all four polymorphs was identical, with a value of about 100 GPa.
Most quantum mechanics programs try to change the length or angle of a chemical bond to simulate a change in structure, since the energy associated with these changes is large.However, the Si-O bond is very mobile and can rotate easily around relatively rigid chemical bonds with very small changes in energy.Thus, changes in crystal structure or dimensional changes in silica are dominated by rotation of torsional bond angles that are very difficult to reproduce with programs such as CASTEP, unless a crystal structure is specified in advance to direct the calculations.
The failure of molecular mechanics or quantum mechanics to predict the volumetric properties of silica polymorphs requires a different modelling technique that can predict properties associated with the relatively small changes in energy associated with physical (Van der Waals) bonding between atoms that are not chemically bonded in a structure.This type of bonding dominates the properties of polymers, since the preferred mode of deformation is through the weak van der Waals bonding normal to the chemical bonding in the chain axis of the polymer macromolecules.A direct analogy is made between a polymeric chain and linked chains of silica tetrahedra.

Potential function method
A modelling technique has been developed to predict the structural properties of polymers, called Group Interaction Modelling, GIM [2].This technique has been used to predict equations of state for polymer-based materials [3].The technique uses an empirical equation to describe the relationship between energy and separation distance between adjacent (but not chemically bonded) groups of atoms in a molecular structure.Fortunately, a very simple power-law function called the Lennard-Jones function works very well for most polymers that are based upon carbon or silicon lattices.To a first approximation, the method assumes that chemical bonds do not deform significantly, relative to the weak physical bonds, due to electronic interactions between atoms or molecules.Note that, although the potential function here is an empirical relation for simplicity, all the fundamental contributions to bond energy are embodied in that function, and a potential function can be constructed from "first principles" by a series of quantum mechanics simulations.
Figure 2 shows the general form of a potential energy well for Van der Waal's bonding.Energy, E, is expressed relative to the depth of the potential energy well, Eo, and dimensions are given as volume, V, relative to the volume at the absolute depth of the potential energy well, Vo.It should be noted that the energies for bonding are negative.
As positive energy (such as thermal energy, HT) is increased, the interaction dimensions move away from the minimum energy position to two possible values on either side of the minimum.Generally, the potential well is asymmetric, such that the equilibrium position of the new well minimum moves to higher volumes to cause thermal expansion in a material.The full method is described in [4] and the pressure can be expressed as a function of volume through the highly non-linear relation where E T and V T are energy and volume at constant temperature T.

Application of model to geological materials
On close inspection, the parameters show a qualitative trend that the reference parameter of volume for the different silica forms scales with the inverse of the density of the polymorphs and the energy parameter stays remarkably constant.The conclusion from this observation is that the many polymorphs of silica can be modelled as a single material in terms of a single adjustable parameter of volume at zero pressure, V T , and two absolute reference parameters of energy and molar volume of an underlying polymeric form of silica, Er and Vr respectively.The key point is that each polymorph adopts a crystal structure as if it were constrained by an energy or pressure, Pr, of configuration associated with the non-zero energy of each particular set of metastable chemical bond angles in the polymorph structure.The remarkably simple conclusion from the pressurevolume relations for silica can be translated into a single predictive equation for any given silica polymorph of the following form This implies that the only variable required to predict the equation of state (EOS) is the initial density of the material.The prediction using this approach for the pressure/volume relationship for each of the crystal forms compared to experimental data is shown in Fig. 3.The results are highly encouraging which gives great confidence in the approach.The model has also been used to predict shock Hugoniot data for a range of geological materials as shown in Fig. 4, again with very encouraging results.

Development of model for sand
The sand used in this study is based on sand to a specification of EN13139 with a particle size distribution (PSD) as shown in Fig. 6 and angular shaped particles.
The Porter-Gould model is capable of predicting the shock EOS response of sand which is dominated by an elastic molecular response at higher impact velocities, as shown in Fig. 7.However, at lower impact velocities the model over predicts the behaviour.
The reason for this is that at the lower impact velocities the behaviour of the sand is dominated by compaction, which is a micromechanical response based around reorganisation of sand grains, fracture and potentially grain plasticity.The compaction response is different for different sands due to particle shape and PSD.At present there is no reliable theoretical method for predicting the compaction response given an initial PSD and shape.Therefore some simplifying assumptions are required.
As a first approximation an analytic approach was taken where it was assumed that all of the complexity of compaction for a specific sand could be represented by a single parameter.The idea being to have a highrate, high pressure parameter that was characteristic of the sand in question and to avoid having to produce a model that required significant number of characterisation experiments.Equation (3) shows how the parameter was calculated.The compaction parameter, C, determines how much of the density increment is inelastic and how much is elastic: that is how much the reference density changes in the transition to the shock state and therefore what the pressure in the elastic equation of state should be.
It was decided that a single plate impact experiment could be used to measure C and then use that value to predict all of the other shock experiments.This is shown in Fig. 8.It can be seen that this analytic approach gives good agreement with the whole data set.Whilst this approach succeeded analytically and showed that there was a potential way forward should C be able to be calculated knowing the microstructure, plastic behaviour and friction coefficients for the sand there was a problem with implementing it into the hydrocodes for which it was meant.The available implementation  schemes could not be structured to use this finite deformation approach and so another way to develop the model for numerical simulations was needed.
The approach adopted for the hydrocode model is to fit the compaction behaviour directly to Quasi-Static data using the Multi-Axial Compression of Concrete at Elevated Temperature (MACET) triaxial rig at Sheffield University.The MACET rig, Fig. 9, comprises six platens which can be controlled independently in both load and displacement control which gives access to many triaxial load paths.
This approach has the advantage of using the current models with updated physical understanding whilst still only needing low rate data to access the high-rate high pressure regime.The EOS QS compaction curve fitted to the MACET data is shown in Fig. 10 for an initial density of 1.54 g/cc.The curve from the MACET data is extended to lock-up on the quartz Hugoniot which has been predicted using the Porter-Gould model.The two points chosen in the region not covered by data are based upon a monotonically increasing bulk modulus.The initial point after yield is chosen so that any realistic starting sand density can be used.
The constitutive response is based on a Johnson-Holmquist model [5] and fitted to QS data as shown in Fig. 11 for dry sand.

Validation of model
The model was implemented in the Eulerian hydrocode GRIM and the Lagrangian hydrocode DYNA.The validation of the model was conducted against experiments in different stress states and regimes.The validation experiments were the plate impact, Split Hopkinson Pressure Bar (SHPB) and fragment penetration.

Plate impact
This is a standard test for measuring the shock response of the material and the experimental arrangement is shown in Fig. 12 and the experiments were performed on the single stage gas gun at the Cavendish Laboratory at Cambridge University.
This particular configuration is also used to measure the release behaviour of the sand under shock loading.The response for dry sand is shown in Fig. 13 and the release path follows the slope of the solid quartz Hugoniot, which is also similar to the assumptions in the model.
Given that there are various assumptions in generating these curves it was felt more realistic to compare the model with the raw heterodyne velocimetry (hetv) laser interferometer data.This comparison is shown in Fig. 14  The level of agreement is encouraging.The arrival time is very close indicating that the wavespeed is correct and also the general structure of the trace is reproduced, particularly the first peak.This suggests that the assumptions in the model concerning the release behaviour are reasonable.

Fragment penetration
The experiments comprise various fragments such as sphere, cubes and Fragment Simulating Projectiles (FSP) impacting finite thickness sand targets (i.e.50 to 200 mm) at a range of velocities (i.e. up to 1500 m/s).The penetration process is illustrated for two impact velocities for a steel cube as shown in Fig. 15.
At higher velocities the fragment is eroded which affects the presented area during penetration and hence the residual velocity.The comparison of the model with experiment for a range of impacts is shown in Fig. 16, which shows excellent agreement.
For semi-infinite penetration the model is within the experimental spread for a range of impact conditions and sand moisture contents.

Split Hopkinson Pressure Bar (SHPB)
The SHPB tests were performed at Sheffield University and the sample was confined in a steel jacket such that the deformation was essentially uniaxial strain as opposed to uniaxial stress.The validation is a direct comparison with    The results are encouraging in terms of general stress levels and also suggests that the wavespeed in the material is correct.

Discussion
The theme of this paper was to develop a model for a new soil when there is a lack of experimental data to fit.
As discussed the major obstacle to developing a model is the compaction behaviour, given that the elastic response is well predicted by the Porter-Gould density-dependent equation of state.This compaction is a function of particle size and particle shape and there are significant barriers to predicting this behaviour including physical properties of the sand, understanding the evolution of sand grain re-arrangement under hydrostatic and other loads, and implementation of relevant schemes into hydrocodes.
The previous sections have demonstrated that the existing models fitted to a QS compaction curve work really well over a range of regimes including shock response and high strain rate behaviour.
The idea therefore is to generate a family of compaction curves for a range of sands and soil 04024-p.5 EPJ Web of Conferences incorporating different particle sizes and shapes.In addition this would need to include moisture content.
Thus for the new material one would use this data to provide the closest fit for the compaction data.This would allow a rapid assessment to be performed for the new material.One issue is that it is not known how many different compaction curves are required to capture the spectrum of response for granular materials.
In terms of developing better theoretical tools the Porter-Gould model is being used to assess compaction behaviour.It has been shown that the elastic response is well predicted by the density-dependent equation of state and so, given experimental data on the response of sand along various load paths it should be possible to deconvolute the inelastic components of the response and thereby quantify both the elastic (i.e.reversible) and inelastic (i.e.irreversible) components.Primary amongst this needed set of data is true hydrostatic compaction of sand coupled to X-ray to understand the trade-off between grain rearrangement, grain fracture and grain plasticity in this highly constrained loading mode.This can then be extended to other interesting loading modes such as isochoric shear as another way in to understanding how these models may be used in current hydrocodes that explicity split hydrostats from deviators.If this, then, allows 1-D deformation to be predicted all the better.
Despite the generally small strains found in geological materials it is considered best to implement this model using a finite deformation scheme which would allow tracking of the inelastic components and leave the stress response to be a hyperelastic outcome.

Conclusions
A sand model has been created based on a theoretical Porter-Gould approach for the EOS combined with a Johnson-Holmquist constitutive model.
The biggest issue is compaction which is a strong function of particle size distribution and shape.
The compaction is fitted to a QS triaxial test (i.e.MACET).
The model has then been shown to predict the behaviour of controlled sand over a wide range of regimes ranging from shock to penetration behaviour giving confidence that this approach works.
It is recommended to generate a family of compaction curves together with additional theoretical developments on separating the reversible and irreversible compaction effects.

Figure 1 .
Figure 1.Different Crystal Forms of Silica.

Figure 2 .
Figure 2. A Potential Well for Intermolecular Energy.

Figure 3 .
Figure 3.Comparison of Model with Experiment for Different Crystal Forms of Silica.

Figure 4 .
Figure 4. Comparison of Model with Shock Hugoniot Data.

Figure 5 .
Figure 5.Comparison of Model with Elastic Unloading Paths.

Figure 7 .
Figure 7. Shock Response Data for Various Dry Sands.

Figure 8 .
Figure 8.Comparison of Compaction Model with Experiment.

Figure 13 .
Figure 13.Shock and Release Paths in Dry Sand.

Figure 14 .
Figure 14.Comparison of Model and Experiment for Het-v Trace.

Figure 16 .
Figure 16.Comparison of Residual Velocity for Fragment.

Figure 17 .
Figure 17.Comparison of Model with Output Trace in SHPB. .