The effect of degree of saturation of sand on detonation phenomena associated with shallow-buried and ground-laid mines

A new materials model for sand has been developed in order to include the effects of the degree of saturation and the deformation rate on the constitutive response of this material. The model is an extension of the original compaction materials model for sand in which these effects were neglected. The new materials model for sand is next used, within a non-linear-dynamics transient computational analysis, to study various phenomena associated with the explosion of shallow-buried and ground-laid mines. The computational results are compared with the corresponding experimental results obtained through the use of an instrumented horizontal mine-impulse pendulum, pressure transducers buried in sand and a post-detonation metrological study of the sand craters. The results obtained suggest that the modified compaction model for sand captures the essential features of the dynamic behavior of sand and accounts reasonably well for a variety of the experimental findings related to the detonation of shallow-buried or ground-laid mines.


Introduction
Detonation of the shallow-buried and ground-laid mines and the subsequent interactions of the resulting shock waves, detonation products and the soil ejecta with the surrounding media and structures involve numerous highly non-linear phenomena of a transient nature.In order to maximize the destructive effects of the explosion or to devise means/strategies for minimizing such effects, a large range of diverse physical phenomena must be considered.While, in principle, one would prefer to study the aforementioned detonation phenomena using an analytical technique, in hope of elucidating the underlying physics of the problem, analytical methods typically entail major simplifying assumptions so that their predictions are often questionable or even contradicted by the experimental observations.Consequently, a better understanding of the explosion phenomena is being gradually gained by combining physical experiments with numerical modeling techniques.This approach is utilized in the present work in which, for example, the experimental results associated with the explosion of shallow-buried and ground-laid C4 mines obtained through the use of an instrumented horizontal mine-impulse pendulum reported in Ref. [1] are compared with a detailed numerical modeling of the same physical problem using AUTODYN, a general purpose non-linear dynamics simulation software [2].
In our recent work [3], a detailed comparison was made between the experimental results reported in Ref. [4] and their computational counterparts for a number of detonation-related phenomena such as the temporal evolutions of the shape and size of the over-burden sand bubbles and of the detonation-products gas bubbles, the temporal evolutions of the side-on pressures in the sand and in air, etc.It was found that the most critical factor hampering a better agreement between the experiment and computational analysis is an inadequacy of the current materials models for sand to capture the dynamic response of this material under blast loading conditions.Hence, the main objective of the present work is to improve the compaction materials model [5] for sand currently implemented in AUTODYN.Specifically, the new materials model for sand is developed to include the effects of saturation and the deformation rate on the constitutive response of the sand.A review of the literature reveals that the currently available materials model for sand suffers from either an inability to account for the rate-dependence of the material's response in the case of the saturated sand [5] or contain a large number of parameters whose estimation via analytical/numerical analysis and /or experimental measurements is very difficult and cumbersome [6][7][8][9][10].
The organization of the paper is as follows.A brief overview of the design, construction and utilization of the instrumented horizontal mine-impulse pendulum is given in Section 2.1.The non-linear dynamics approach, the relevant materials models and the definition of the computational problem investigated are respectively discussed in Sections 2.2, 2.3 and 2.4.The results obtained in the present work are presented and discussed in Section 3. The main conclusions resulting from the present work are summarized in Section 4.

An overview of the horizontal mine-impulse endulum experiment
Since one of the objectives of the present work is to assess the validity of the newly proposed materials model for sand by comparing the present non-linear dynamics based computational results with their experimental counterparts obtained using an instrumented horizontal mine-impulse pendulum used in Ref. [1], a brief overview of the construction and experimental procedure associated with the impulse pendulum is presented in this section.A more detailed account of the design and the construction of the instrumented horizontal mine-impulse pendulum can be found in Ref. [11].
The instrumented horizontal mine-impulse pendulum consists of a 5 m long horizontal steel arm with a 1200 mm × 1200 mm square measuring pan placed at the free end of the arm 400 mm above the ground.The arm is attached to the base assembly at the other end through a horizontal pivot.The charge, typically consisting of a cylindrically shaped (14.6 cm in diameter and 5 cm high) C4 mine is placed under the center of the measuring pan and detonated.The mine is either laid on the ground or is buried to different depths.The resultant maximum angular displacement of the pendulum arm is measured and used to calculate the detonation-induced impulse on the pendulum, see Eq. ( 5).The use of the mine-impulse pendulum enabled an investigation of the effects of the sand type/properties, extent of saturation with water, the target stand-off distance and the mine depth of burial on the total detonation-induced impulse [1].In an earlier design, the measuring pan was constructed of mild steel, however, the initial experiments revealed that such a measuring pan undergoes substantial plastic deformation.Consequently, the central 600 mm × 600 mm section of the measuring pan was replaced with a 50 mm thick Rolled Homogenized Armor (RHA) plate.The deformation of the RHA plate was found to be in the order of 10 −4 to 10 −3 which is small and thus justifies the assumption of negligible deformation.The maximum angular deflection of the pendulum was obtained using a combination of the following three methods: (a) a cable potentiometer, (b) a mechanical gage and (c) high speed video recording of a large pointer.
The relationship between the total detonation-induced impulse on the pendulum and the maximum angular displacement of the pendulum mentioned above was derived in Ref. [11] under the assumption that both the effect of gravity and that of pendulum displacement during the initial loading phase (the time period over which the momentum is transferred from the detonation-products/sand ejecta to the measuring pan) can be neglected.These were justified by the fact that a typical duration of the initial loading phase is about 1-2 ms, while the typical time for the pendulum to reach the maximum angular displacement position is around 300-400 ms.A brief overview of the derivation of the impulse vs. maximum angular displacement relation is given below.The derivation of this relation is based on the principle of conservation of the total energy (a sum of the potential and kinetic energies).When applied to the initial position of the pendulum and the position of the pendulum corresponding to the maximum angular displacement, this principle leads to the following equation: where I 0 is the moment of inertia of the pendulum arm, ω 0 is the initial angular velocity of the pendulum arm, m is the mass of the pendulum arm, g is the gravitational acceleration (9.81 m/s 2 ), r is the distance between the pivot point and center of gravity of the arm and θ max is the maximum angular displacement of the pendulum arm.It should be noted that the potential energy of the pendulum in its initial position is arbitrarily set to zero (the left hand side of Eq. ( 1)), while the kinetic energy of the pendulum associated with the maximum angular displacement is also zero (the right hand side of Eq. ( 1)).Thus the left hand side of Eq. ( 1) defines the initial rotational kinetic energy of the pendulum while the right hand side of the same equation defines the maximum potential energy of the pendulum.
It should be noted that it is postulated in Eq. ( 1) that the kinetic energy is initially imparted to the pendulum by the detonation products/soil ejecta without any pendulum movement.
The initial angular velocity of the pendulum, ω 0 , is obtained from the following angular impulse relationship: where F is the detonation-induced normal force acting on the pendulum over a time period between t 0 and t f , t is the time and R is the distance between the point of application of the force and the pivot point.Since the blast loads act normal to the surface of measuring pan and over a very short time period, R can be considered as a constant.
If the effect of gravity is neglected on the right hand side of Eq. ( 2) and a use is made of the definition of the total detonation-induced impulse on the pendulum, J: the initial angular velocity can be defined as: Substitution of Eq. (4) in Eq. (1) yields Thus for a given pendulum with the parameters m (= 1480 kg), I 0 (= 14,700 kg-m 2 ), r (= 2.44 m) and R (= 4.27 m) Eq. ( 5) enables determination of the detonation-induced impulse from the measured values of the maximum angular displacement (θ max ) of the pendulum.The values for the pendulum parameters listed above correspond to the ones for the horizontal mine-impulse pendulum used in Ref. [1].

Non-linear dynamics modeling of detonation phenomena
All the calculations carried out in the present work are done using AUTODYN, a general purpose non-linear dynamics modeling and simulation software [2].In this section, a brief overview is given of the basic features of AUTODYN, emphasizing the aspects of this computer program which pertain to the problem at hand.
AUTODYN is a fully integrated engineering analysis computer code which is particularly suited for modeling explosion, blast, impact and penetration events.Codes such as AUTODYN are commonly referred to as "hydrocodes".Within the code, the appropriate mass, momentum and energy conservation equations coupled with the materials modeling equations and subjected to the appropriate initial and boundary conditions are solved.The numerical methods used for the solution of these equations involve finite difference, finite volume and finite element methods and the choice of the method used (i.e."processor" as referred to in AUTODYN) depends on the physical nature of the problem being studied.The power of AUTODYN is derived mainly from its ability to handle complex problems in which different regions can be analyzed using different methods such as the Lagrange processor (typically used for solid continuum and structures) and the Euler processor (commonly used for modeling gases, liquids or solids subject to large deformations).While the available Euler processor provides multi-material capabilities, an additional Euler-FCT single material processor in which materials are combined to a single material using a Flux Corrected Transport (FCT) approach is available to help handle computationally intensive multi-material blast phenomena.
Additional methods available in AUTODYN include: an ALE (Arbitrary Lagrange Euler) processor capable of carrying out an automatic rezoning of distorted grids; a Shell processor designated for modeling thin structures and a gridless SPH (Smooth Particle Hydrodynamics) processor which does not suffer from a grid tangling problem (typically encountered in Lagrangian processor) and does not entail the use of an unphysical erosion algorithm (removal of highly distorted grids to help the numerical procedure).The authors are not aware of any comprehensive open literature on the study of the effect of solver type on the computational results in problems which can be solved with alternative solvers.AUTODYN reference manual provides a comparison between the results obtained using Euler, Lagrange and ALE solvers for the case of a cylindrical projectile impacting a thick plate and shows that the solver choice has a relatively small effect on the computational results.
In the present work, the Euler-FCT processor was used to represent air and C4 gaseous detonation products.Sand is handled using the Lagrange processor while the various components of the instrumented horizontal mine-impulse pendulum are represented using a shell processor.The interactions between the different processors are accounted for through the use of the part-interaction options within AUTODYN [2].

Materials constitutive models
Hydrodynamic computer programs such as AUTODYN [2] are capable of predicting an unsteady, dynamic motion of a material system by solving the appropriate mass, momentum and energy conservation equations, subjected to the associated initial and boundary conditions.However, for the aforementioned boundary value problem to be fully defined, additional relations between the flow variables (pressure, density, energy, temperature, etc.) have to be defined.These additional relations typically involve an equation of state, a strength equation and a failure equation for each constituent material.These equations arise from the fact that, in general, the total stress tensor can be decomposed into a sum of a hydrostatic stress (pressure) tensor (which causes a change in the volume/density of the material) and a deviatoric stress tensor (which is responsible for the shape change of the material).An equation of state then is used to define the corresponding functional relationship between pressure, density and internal energy (temperature), while a strength relation is used to define the appropriate equivalent plastic-strain, equivalent plastic strain rate, and temperature dependences of the equivalent deviatoric stress.In addition, a materials model generally includes a failure criterion, i.e. an equation describing the (hydrostatic or deviatoric) stress and/or strain condition which, when attained, causes the material to fracture and loose its ability to support normal and shear stresses.
In the present work the following four materials are utilized within the computational domain: air, sand, AISI 1006 steel and Rolled Homogenized Armor (RHA).In the following sections, a brief description is given of the models used for each of the four constituent materials.The values of all the material parameters defined in the remainder of the section are available in the AUTODYN materials library [2].The data cannot be disclosed here due to copyright violation concerns.

Air
Air is modeled as an ideal gas and, consequently, its equation of state is defined by the ideal-gas gamma-law relation as [2]: where P is the pressure, γ the constant-pressure to constant-volume specific heats ratio (= 1.4 for a diatomic gas like air), ρ 0 (= 1.225 kg/m 3 ) is the initial air density, and ρ is the current density.For Eq. ( 6) to yield the standard atmosphere pressure of 101.3 kPa, the initial internal energy E is set to 253.4 kJ/m 3 which corresponds to the air mass specific heat of 717.6 J/kg•K and a reference temperature of 288.2 K.
Due to the use of a single-material Euler-FCT processor for the gas-phase region, the C4 detonation products are not modeled as a separate material within the gas phase.Rather, C4 detonation products are modeled initially as a cylindrically shaped air region with a high density ρ (= 1601 kg/m 3 ), and a high internal energy density, e (= 5.621•10 6 J/kg).The corresponding detonation products pressure and the fire ball temperature take on values of 18 MPa and 2950 K, respectively.
Since air is a gaseous material and has no ability to support either shear stresses or negative pressures, no strength or failure relations are required for this material.

AISI 1006 steel
With the exception of the 600 mm × 600 mm central square section of the measuring pan, the instrumented horizontal mine-impulse pendulum was constructed from medium-carbon AISI 1006 Steel.For inert solid materials like AISI 1006 steel a linear type of equation of state is typically used which assumed a Hooke's law type relationship between the pressure and the volume change as: where K is the bulk modulus of the material and µ = ρ ρ0 − 1 is the compression ratio.Within AUTODYN material database, the initial material density ρ 0 , the bulk modulus K, the specific heat and the reference temperature are defined for AISI 1006 steel.
To represent the constitutive response of AISI 1006 steel under deviatoric stress, the Johnson-Cook model is used.This model is capable of representing the material behavior displayed under large-strain, high deformation rate, high-temperature conditions, of the type encountered in problems dealing with the interactions of detonation products and solid structures.Within the Johnson-Cook model, the yield stress is defined as: where ε pl is the equivalent plastic strain, εpl is the equivalent plastic strain rate, A is the zero plastic strain, unit plastic strain rate, room temperature yield stress, B is the strain hardening constant, n is the strain hardening exponent, C 1 is the strain rate constant, m is the thermal softening exponent and T H = (T − T room )/(T melt − T room ) is a room temperature (T room ) based homologous temperature while T melt is the melting temperature.All temperatures are given in degrees of Kelvin.Since the sections of the instrumented horizontal mine-impulse pendulum constructed from the AISI 1006 steel are generally subjected to relatively low stresses, no failure model was used for the AISI 1006 steel in the present work.

Rolled Homogenized Armor (RHA)
As mentioned earlier, the 600 mm × 600 mm central square section of the measuring pan was constructed from a Rolled Homogenized Armor (RHA) plate material.The same type of materials models (a linear equation of state and a Johnson-Cook strength model) used for the AISI 1006 steel are also used to represent the dynamic response of the RHA plate material.However, the values of the model parameters differ for the two typed of materials.

Sand
Sand has generally a complex structure consisting of mineral solid particles which form a skeleton.The pores between the solid particles are either filled with effectively dry air (this type of sand is generally referred to as "dry sand"), with water ("saturated sand") or a two-phase water/air mixture ("unsaturated sand").The relative volume fractions of the three constituent materials in the sand (the solid mineral particles, water and air) are generally quantified by the porosity, α, and the degree of saturation (Saturation Ratio), SR, which are respectively defined as and where V p is the volume of void (pores), V w is the volume of water and V is the total volume.
Surface roughness and the presence of inorganic/organic binders are generally considered to be the main causes for friction/adhesion at the inter-particle contacting surfaces.Deformation of the sand is generally believed to involve two main basic mechanisms [12]: (a) elastic deformations (at low pressure levels) and fracture (at high pressure levels) of the inter-particle bonds and (b) elastic and plastic deformations of the three constituent materials in the sand.The relative contributions of these two deformation mechanisms as well as their behavior are affected primarily by the degree of saturation of sand and the deformation rate.Specifically, in dry sand the first mechanism controls the sand deformation at low pressures while the second mechanism is dominant at high pressures and the effect of deformation rate is of a second order.In sharp contrast, in saturated sand very low inter-particle friction diminishes the role of the first deformation mechanism.On the other hand, the rate of deformation plays an important role.At low deformation rates, the water/air residing in the sand pores is squeezed out during deformation and, consequently, the deformation of the sand is controlled by the deformation of the solid mineral particles.At high pressures, on the other hand, water/air is trapped within the sand pores and the deformation of the sand is controlled by the deformation and the volume fractions of each of the three constituent phases.
Within AUTODYN, the dynamic response of sand is represented using a compaction materials model which was formulated using the experimental results obtained by Laine and Sandvik [5].A brief description of the compaction materials model is given below.

Compaction Materials Model for Sand
The "compaction" equation of state for sand is based on a piece-wise linear pressure-density relation schematically shown in Fig. 1.It should be noted that, since pressure does not depend explicitly on the internal energy, this relation is equivalent to the standard Mie-Gruneisen equation of state in which the Gruneisen gamma parameter, Γ = v ∂P ∂E v is set to zero.This means that the model would give a more reliable material response under the conditions when either the energy absorbed is not very high (e.g. when the applied pressure levels are not significantly larger than the pressure levels at which the porous material crushes and compacts into a solid material), when the initial material porosity is small or when the magnitude of the Gruneisen gamma parameter is near zero.
Within the AUTODYN computer program [2], the initial density of the porous material at a zero pressure level, P 1 is denoted as ρ 1 .As the pressure is applied, the relation between the pressure and the density (denoted in Fig. 1 as "Plastic Compaction") is defined using up to ten (ρ, P ) pairs of values.This portion of the pressure vs. density relation is associated with a permanent, plastic compaction of sand.Full compaction of the sand corresponds to the last pair of the (ρ, P ) values of the plastic compaction curve.
An increase in pressure beyond the point of full compaction is defined by the following elastic loading linear pressure-density relation: where C s is the sound speed in fully compacted sand at zero pressure and ρ s is the mass density of the fully-compacted sand under a zero applied pressure.Elastic unloading/reloading of the porous material like sand at any level of compaction is generally governed by the following differential equation: where C is the sound speed in sand at a density ρ.As indicated in Fig. 1 by the curves denoted as "Elastic Unloading/Reloading", the pressure-density relation during elastic unloading/reloading is not linear which is due to the fact that the sound speed in sand is a function of the material density.Within AUTODYN [2], density dependence of the sound speed, C, is defined as a piece-wise linear relation in terms of up to ten (ρ, C) pairs of values.
The "compaction" strength model for sand is based on an isotropic, perfectly plastic, rate independent yieldsurface approximation and postulates that the yield stress depends explicitly on pressure and not on material density.Within the AUTODYN program [2], the relationship between the yield stress, Y , and pressure, P , is defined as a piece-wise linear function consisting of up to ten (P, Y ) pairs of values.The yield stress quantifies the resistance of the material to a plastic (irreversible) shape change.The plastic shape change occurs when the magnitude of the second invariant of the deviatoric part of the stress tensor becomes equal to the yield stress.
Unloading (and subsequent reloading) of a previously plastically deformed material is of an elastic (reversible) nature and, in this case, the deviatoric stress is proportional to the deviatoric strain with the proportionality constant being equal to the shear modulus, G.In a porous material such as sand, the shear modulus is a function of the material density.The "compaction" G vs ρ. relation is defined within AUTODYN [2] as a piece-wise linear function using up to ten (ρ, G) pairs of data.
The failure behavior of sand is modeled within the AUTODYN materials database by specifying a minimum (negative) value of the hydrodynamic pressure below which, the material fractures, and looses its ability to support any tensile or shear stress.However, if a given "fractured" material region is subsequently subjected to positive pressures, it is given an ability to reheal and close up its cracks.

Modified Compaction Materials Model for Sand
The compaction materials model for sand described in the previous section does not include two very important factors controlling the dynamic response of sand under the blast loading conditions i.e., the effects of the degree of saturation and the deformation rate.In this section, an attempt is made to modify the original compaction materials model for sand in order to incorporate these two effects.The modifications of the compaction materials model for sand, proposed in the present work, are based on the following set of assumptions: (a) The average sand particle size, particle size distribution and the presence of inorganic/organic natural matter in sand all have a second order effect on the dynamic constituent response of the sand.This assumption was justified by the experimental observations reported in Ref. [1] which clearly showed that the effect of sand type (e.g.prairie sand containing high level of silt and clay, impurity-free −30/+50 sand, etc.), on the detonation-induced momentum transfer to the instrumented horizontal mine-impulse pendulum was small in comparison with the effect of the degree of saturation; (b) The dynamic mechanical response of the sand at any degree of saturation can be obtained as a linear combination of the corresponding dynamic material behaviors for the dry and the saturated sand; (c) The dynamic mechanical response of the dry sand is not rate dependent and it can be represented by the original compaction model implemented in AUTODYN [2]; (d) The initial density of the saturated sand, ρ 1,sat , can be calculated using densities of the solid material in the sand, ρ s , and water, ρ w , and the known level of sand porosity, α = 1 − (ρ 1 /ρ s ), as: (e) When the saturated sand is subjected to relatively low deformation (compression) rates, water is given enough time to leave the pores and hence the density of the fully compacted sand and the pressure at which full compaction is attained are identical to their counterparts in the dry sand; (f) When the saturated sand is subjected to very high compression rates, water will be trapped inside the pores and, due to a very low compressibility of the water, the compressibility of the sand is controlled by the compressibility of its solid phase.In other words, the saturated sand behaves as a fully compacted sand under high deformation rates and can only undergo an elastic compaction; (g) Under intermediate deformation rates, the dynamic material response of the saturated sand can be obtained using a linear interpolation of the high-low-deformation rate behaviors of the saturated sand.A value of 1.0•10 5 s −1 is used as the "high" deformation rate, εhigh , and a value of 1.0•10 −3 s −1 is used as the "low" deformation rate, εlow .At the deformation rates exceeding 1.0•10 5 s −1 and at the deformation rates below 1.0•10 −3 s −1 , the dynamic behavior of sand is assumed to be rate independent and to correspond to the dynamic sand behavior at the respective (1.0•10 5 s −1 or 1.0•10 −3 s −1 ) deformation rates.The linear interpolation of the dynamic sand behavior at the intermediate deformation rates was based on the logarithms of the deformation rates as: where the densities ρ, ρ high and ρ low correspond respectively to the deformation rates ε, ε high and εlow and are all associated with the same level of pressure.The computational results obtained are found not to be significantly affected by an order of magnitude changes in the values for the high and low deformation rates; and (h) Since the irreversible plastic deformation of sand is dominated by the plastic deformation behavior of its solid phase, it is assumed to be independent of the degree of saturation and the rate of deformation.In other words, only the equation of state in the original compaction model was modified following the aforementioned procedure.The pressure vs. density relations corresponding to the original compaction model for sand and also for dry sand at a porosity level of 38% is displayed in Fig. 1.The pressure vs. density relations corresponding to saturated sand in the present formulation with a porosity level of 38% is displayed in Fig. 2.
The new equation of state is implemented in the user subroutine "mdeosuser.f90" and interfaced with the AUTODYN computer program [2].

Problem definition and computational analysis
In this section, a brief description is given of the computational model used to simulate the interaction of the detonation-products/soil ejecta resulting from the explosion of a shallow-buried or ground-laid mine and the instrumented horizontal mine-impulse pendulum.The computational modeling of this interaction involved two distinct steps: (a) geometrical modeling of the instrumented horizontal mine-impulse pendulum and (b) a non-linear dynamics analysis of the momentum transfer from the detonation-products/soil ejecta to the pendulum.
Various computational domains used in the present study are shown in Fig. 3.The geometrical models for the various components of the pendulum were constructed using 50 mm × 50 mm square shell elements.An advantage was taken of the planar symmetry of the model.In other words, a vertical plane of symmetry was placed along the length of the pendulum which enabled only a half of the pendulum to be modeled.In accordance with the instrumented horizontal mine-impulse pendulum used in Ref. [1], different sections of the pendulum were constructed using AISI 1006 steel and RHA plate material.Welded joints of the different sections of the pendulum were simulated by joining the components in question.
The head of the pendulum was placed in an Euler-FCT region consisting of 74,000 25 mm edge-length cubic cells.In the case of a surface laid mine, the mine was represented by a high-density high-energy cylindrical air region located within the Euler-FCT domain.In the case of a shallow-buried mine, two joined Lagrange domains were used to define a sand region containing a cylindrical cavity whose shape and size match those of the C4 mine (diameter 132 mm and height 35 mm).A second Euler-FCT domain overlapping with the two sand domains is defined and the portion of this domain corresponding to the cylindrical sand cavity defined above is initially filled with high-density high-energy air.The Euler-FCT is capable of capturing the shock development in air.
The air/sand and air/pendulum interactions are accounted for using the appropriate Euler/Lagrange coupling option with AUTODYN [2].Likewise, the sand/pendulum interactions were modeled through the use of the appropriate Lagrange/Lagrange coupling option.
At the beginning of the simulation, the pendulum is assumed to be at rest (with the gravitational force acting downwards), while the Lagrange and Euler-FCT domains are filled with stationary materials (sand and air, respectively).As explained earlier, the C4 mine was initially modeled as a cylindrical high-density, high-energy sub-domain within the Euler-FCT region.
The motion of the pendulum was constrained to within a vertical plane and a fixed single-point constraint was applied to its pivot point.The "flow out" boundary conditions were applied to all the free faces (the faces which do not represent interfaces between the different domains) of the Euler-FCT domain except for the face associated with the vertical symmetry plane.To reduce the effect of reflection of the shock waves at the outer surfaces of the Lagrange domain, "transmit" boundary conditions were applied to all the free faces of this domain except for the face associated with the vertical symmetry plane.
To speed up the calculations, all Euler-FCT and Lagrange domains were removed from the analysis after approximately 10 ms following detonation when the extent of interaction between the detonation-products/sand ejecta and the pendulum was negligibly small.No numerical difficulties were encountered during implementation of this procedure and the computed results were found to be essentially identical when the Euler-FCT and Lagrange domains were removed at a simulation time of 50 ms.

Implementation validation of the modified compaction equation of state
To validate that the modified compaction equation of state, as described in Section 2.3.4, is correctly implemented in the mdeosuser.f90user subroutine and correctly interfaced with AUTODYN [2], a series of one cubic-element hydrostatic-compression analyses under different deformation rates is carried out in the present section and the results of these analyses displayed using pressure vs. density plots.The corresponding plots were also generated via numerical integration of the modified compaction equation of state using MATLAB [13], a general purpose mathematical package.A perfect agreement was found between the two sets of pressure vs. density plots confirming that the implementation of the modified compaction equation of state was correct.
An example of the results obtained during the validation of the implementation of the equation of state are shown in Fig. 4(a) and (b).It should be noted that the results obtained using AUTODYN [2] and MATLAB [13] are indistinguishable from each other.

Detonation-induced momentum transfer to the pendulum
The non-linear dynamics analysis described in Section 2.4 is carried out in the present section in conjunction with the new materials model for sand described in Section 2.3.4 and Eq. ( 5) to determine the total momentum (impulse) transferred to the horizontal instrumented mine-impulse pendulum by the detonation-products and the sand ejecta.The results obtained are compared with their experimental counterparts as reported in Ref. [1].
The effect of the degree of saturation in sand on the total impulse transferred to the pendulum for two types of sand (prairie sand containing high level of silt and clay and impurity-free −30/+50 mesh sand) for the case of a shallow-buried mine at a 5 cm depth of burial obtained in Ref. [1] is displayed in Fig. 5. Also displayed in Figure 5 are the computational results obtained in the present work using both the original and the modified materials model for sand.The results displayed in Fig. 5   (a) The type of sand appears to have a relatively small effect on the impulse transfer.This conclusion was partly based on the results displayed on Fig. 5. and partly on the conclusions made in Ref. [1].(b) The degree of saturation of sand has a major effect on the impulse transfer increasing it by 100-150% in saturated sand relative to that in nearly dry sand; (c) The original compaction model neglects the effect of moisture on the materials response of the sand and hence the total impulse transferred to the pendulum is independent of the degree of saturation, which is consistent with the experimental observation; (d) The modified compaction model, on the other hand, predicts an increase in the total impulse transferred to the pendulum with an increase in the degree of saturation of sand in a qualitative agreement with experimental data reported in Ref. [1]; and (e) The quantitative agreement between the computational total impulse vs. degree of saturation results obtained using the modified compaction model and the experimental results is significantly improved with respect to the agreement between the original model and the experimental results.
The effect of the degree of saturation in sand on the total impulse transferred to the pendulum in the case of the prairie sand containing high level of silt and clay for four different locations of the C4 mine obtained in Ref. [1] is displayed in Fig. 6(a) and (b).To help the reader better understand the definition of DOB, a simple sketch of various mine locations is included in Fig. 6(b).The location of the mine is denoted by the corresponding value of the 'Depth of Burial' (DOB).A 0 cm DOB corresponds to the flush-buried mine while a −5 cm DOB corresponds to a ground-laid mine.Also displayed in Fig. 6(a) and (b) are the computational results obtained in the present work using the modified materials model for sand.The results displayed in Figures 6(a) and 6(b) can be summarized as follows: (a) The lowest value of the impulse transferred to the horizontal instrumented mine-impulse pendulum is obtained in the case of a ground-laid mine, Fig. 6(a), since this transfer takes place almost exclusively via the interaction of the gaseous detonation products with the pendulum.This is supported by the fact that the impulse transferred to the pendulum does not show a consistent increase with the increase in the degree of saturation; (b) For a flush-buried mine (0 cm DOB), Fig. 6(a), the detonation induced impulse transfer is increased since, in addition to the detonation products, sand ejecta also interact with the pendulum; (c) The largest impulse transfer occurs in the case of shallow-buried mines (5 cm and 10 cm DOB), Fig. 6(b), where the extent of sand ejection and interaction with the pendulum is the largest; (d) Since the total impulse transferred to the pendulum is somewhat larger for the case of 5 cm DOB than that in the case of 10 cm DOB, it appears that there is an optimum DOB which maximizes the lethal effect of detonation of a shallow-buried mine.This can be rationalized by the fact that as the DOB is increased the effects of detonation become more confined within the soil (the camouflet effect); and (e) The overall quantitative agreement between the computational results and their experimental counterparts at different values of DOB is reasonable considering the fact that: (i) as mentioned earlier, the original compaction model is used to represent the behavior of the dry component of the sand and (ii) a significant disagreement between the computational and the experimental results is seen in the −5 cm DOB case, Fig. 6 where the choice of the materials model for sand is essentially immaterial since only the gaseous detonation products interact with the pendulum.The observed discrepancy between the computed and the experimental results in the case of a −5 cm DOB is a measure of the short comings of AUTODYN to accurately account for free-air blast phenomena.

Blast wave propagation in sand
To further test the validity of the modified compaction model for sand, a non-linear dynamics analysis of the blast wave propagation in sand is carried out and the results obtained compared with their experimental counterparts across the main boundaries without reflection, Ref. [10].Since the charge volume is very small in comparison to the volume of the computational domain, the mesh size was increased as a function of the distance from the charge center.To check the adequacy of the mesh size in simulating the blast effects, a convergence study was conducted, whereby trial calculations were performed and the mesh was refined after each analysis until the difference of the results between two consecutive calculations was deemed small.The final mesh used consisted of 35,250 elements, with the smallest element being 25 mm × 25 mm, and the largest element being 200 × 200 mm.A series of gage points were defined within the sand domain at the same depth as the buried charge to monitor the blast wave propagation.The TNT charge was represented using the Jones-Wilkins-Lee (JWL) equation of state, Ref. [14], while the modified compaction model was used for sand.
The variation of the peak pressure with a scaled distance from the charge center in sand at four different Degrees of Saturation (DOS) obtained experimentally in Ref. [10] is displayed using a log-log plot in Fig. 7(a) and (b).For comparison, the corresponding computational results obtained in the present work are also shown in Fig. 7(a) and (b).The scaled distance is defined as the ratio of the distance from the charge center and the cube root of the charge mass.The results displayed in Fig. 7 The variation of the specific impulse with a scaled distance from the charge center in sand at four different degrees of saturation obtained experimentally in Ref. [10] is displayed in Fig. 8(a) and (b).For comparison, the corresponding computational results obtained in the present work are also shown in Fig. 8(a) and (b).The results displayed in Fig. 8(a) and (b) are obtained by integrating with respect to time the pressure vs. time traces at the location of the pressure transducers/gage points.An analysis of the results displayed in Fig. 8(a) and (b) reveals that conclusions similar to those drawn from Fig. 7(a) and (b) can be made.

Sand crater morphology
The final comparison between the experimental and computational results used to validate the modified compaction model for sand is presented in this section.The experimental results used in this section pertain to the sand crater morphology resulting from the detonation of 0.1 kg C4 charges buried in sand (contained in a thick-wall barrel) at different DOBs obtained in Ref. [4].Since a detailed account of the associated computational procedure used is given in our previous work [3], such details are not presented here.Rather, only the relevant results are presented and compared with their experimental counterparts.
The final morphology of the craters resulting from detonation of the C4 high-energy explosive at 0 cm, 3 cm and 8 cm DOBs experimentally determined in Ref. [4] are displayed in Fig. 9(a)-(c), respectively.Figure 9 (c) For the cases of 0 cm and 3 cm DOB, the central portion of the crater appears to be nearly flat, Fig. 9 The corresponding computational results obtained using the original compaction model for sand are displayed in Fig. 9  A careful examination of the time evolution of the pressure field revealed that formation of the bulged region at the bottom of the crater is a result of the complex interactions between the shock waves and the rarefaction waves in the air above the sand.
The corresponding computational results obtained using the modified compaction model for sand are presented in Fig. 9(g)-(i).A comparison of the results displayed in Fig. 9(a)-(c), (d)-(f) and (g)-(i), show that the modified compaction model for sand gives a qualitatively better agreement with the experimental results than the original compaction model.
It should be noted that the measured results shown in Fig. 9(a)-(c) correspond to the final crater shapes while the computed crater shapes displayed in Fig. 9(d)-(f) and (g)-(i) are obtained after simulation time of 150 ms (at this time, sand velocities are found to be quite small, suggesting that no major subsequent changes in the crater shape should be expected).In other words, the crater shape at an infinite time is assumed to be the same as that obtained at the simulation time of 150 ms.
To obtain a more quantitative comparison between the measured and computed crater shapes, the corresponding variations in the crater depth and the crater width with the charge DOB are displayed in Fig. 10(a) and (b),respectively.It should be noted that the measured crater depths correspond to their final values while the measured crater widths correspond to the time of 12 ms following detonation, the time which was matched in the computational analysis.The crater widths at 12 ms post-detonation time were determined in Ref. [4] using high speed photography.Hence, obtaining a better agreement with respect to the crater width between the experiment and the computational analysis is more critical for the successful validation of the new materials model for sand.
By analyzing Figs 9-10, the following main observations can be made: Pivot point and point of application of normal blast force r -Distance between pivot point and center of gravity of Pendulum arm ρ quantity pl -Plastic state quantity room -Room temperature quantity s -Fully-compacted sand related quantity H -Homologous quantity sat -Saturation related quantity w -Water related quantity Superscripts m -Thermal softening exponent n -Strain hardening exponent

Fig. 1 .
Fig. 1.Pressure vs. density relations for dry sand as defined in the AUTODYN materials library [2].

Fig. 2 .
Fig. 2. Pressure vs. density relations for saturated sand at low and high deformation rates proposed in the present work.

Fig. 3 .
Fig. 3. Various Computational domains used in the present non-linear dynamics analysis.
can be summarized as follows:

Density 1 Fig. 4 .
Fig. 4. Pressure vs. density relation for: (a) fully saturated sand with a porosity level of 38% at various deformation rates and (b) for sand at different degrees of saturation and a deformation rate of 1.10 3 s −1 .

Fig. 5 .
Fig.5.The effect of degree of saturation and the sand type on the total impulse transferred to the instrumented horizontal mine-impulse pendulum.Explosive weight and type: 1.34 kg C4; Depth of burial: 5 cm.

Fig. 6 .
Fig.6.The effect of degree of saturation and the depth of burial on the total impulse transferred to the instrumented horizontal mine-impulse pendulum.Explosive weight and type: 1.34 kg C4; Sand type: Prairie Sand.
(a) and (b) can be summarized as follows: (a) At each level of the degree of saturation, the peak pressure decreases monotonically (non-linearly) with the scaled distance; (b) At any value of the scaled distance, the peak pressure increases with an increase in the degree of saturation of sand; (c) The trends identified in (a) and (b) are also displayed by the computational results; (d) The overall quantitative agreement between the computational and experimental results is reasonable considering a significant scatter in the experimental results.It should be noted that the results displayed in Fig. 7(a) and (b) are a better indication of the ability of the modified compaction model to account for the dynamic behavior of sand at different levels of the degree of saturation than the results presented in Figs 5 and 6 since: (i) they include only the effect of charge/sand interactions and (ii) the detonation products are represented by a more realistic (JWL) model.(e) It is unfortunate that the experimental pressure vs. time traces were not available, since such data would allow for a more definite validation of the new material model for sand.

Fig. 8 .
Fig. 8. Variation of specific impulse with the scaled distance from the charge center and the degree of saturation.Explosive weight and type: 8 kg TNT; Sand type: Sandy loam.
(a)-(b), while for the case of 8 cm DOB, Fig. 9(c), the central portion of the crater contains a minor bulge. (d)-(f).

Figure 9 (Fig. 9 .
Fig. 9. Effect of the Depth of Burial (DOB) on the crater shape obtained: (a)-(c) experimentally in Ref. [4]; (d)-(f) computationally using original compaction model for sand; (g)-(i) computationally using the modified compaction model for sand.Explosive weight and type: 0.1 kg C4; Sand type: −30/+50 mesh silica sand; DOS: 7.5%.(b) The computed results show that some displaced sand remains above the initial position of the sand/air interface (denoted as a dashed line); (c) While the computed results show an increase in the crater depth with an increase in DOB, in agreement with the measured results, this variation is substantially more pronounced in the case of the computed results; (d) The computed values of the crater depth at low values of DOB, Fig. 9(d)-(e), are substantially lower than their measured counterparts, Fig. 9(a)-(b); and (e) While the computed crater shape for the largest DOB, Fig. 9(f), shows a bulge at its bottom in agreement with the corresponding experimentally determined crater shape shown in Fig. 9(c), the height of the computed bulge is clearly smaller.