Discrete Element Simulation of Elastoplastic Shock Wave Propagation in Spherical Particles

Elastoplastic shock wave propagation in a one-dimensional assembly of spherical metal particles is presented by extending wellestablished quasistatic compaction models. The compaction process is modeled by a discrete element method while using elastic and plastic loading, elastic unloading, and adhesion at contacts with typical dynamic loading parameters. Of particular interest is to study the development of the elastoplastic shock wave, its propagation, and reflection during entire loading process. Simulation results yield information on contact behavior, velocity, and deformation of particles during dynamic loading. Effects of shock wave propagation on loading parameters are also discussed. The elastoplastic shock propagation in granular material has many practical applications including the high-velocity compaction of particulate material.


Introduction
The dynamic response of the granular media has become increasingly important in many branches of engineering.It includes material processing involving dynamic compaction and material processing, as well as acoustics and wave propagation in geomechanics.The granular matter shows discrete behavior when subjected to static or dynamic loading [1][2][3].The dynamic wave propagation in granular media shows distinct behavior from the wave propagation in continues media [1].Shukla and Damania [4] discuss the wave velocity in granular matter and shown experimentally that it depends upon elastic properties of the material and on geometric structure.Similarly, Shukla and Zhu [5] investigate explosive loading discs assembly and found that the force propagation through granular media depends on impact duration, arrangement of the discs, and the diameter of discs.Tanaka et al. [6] investigate numerically and experimentally the dynamic behavior of a two-dimensional granular matter subjected to the impact of a spherical projectile.
To investigate dynamic response, many researchers [7][8][9][10][11] have modeled the granular matter as spherical particles using the micromechanical modeling of contact between particles.These studies focused on equivalent macroelastic constitutive constants during dynamic loading.Similarly experimental work using dynamic photoelasticity and strain gage are performed to investigate contact loads between particles both under static and dynamic loading [12,13].Sadd et al. [12] perform numerical simulations to investigate the effects of the contact laws on wave propagation in granular matter.Similarly Sadd et al. [13] use the discrete element method (DEM) to simulate wave propagation in granular materials.Results of this study show wave propagation speed and amplitude attenuation for two-dimensional assembly of spherical particles.However, this study is restricted to the elastic range only while the material stiffness and damping constants used in the model are determined by photoelasticity.DEM was initially developed by Cundall and Strack [14] and this numerical method has been widely used for granular material simulations [15][16][17].Different engineering approaches are discussed in [18,19] to model the behavior of granular matter using DEM.Dynamic compaction of metal powder is also reported in the literature [20][21][22][23][24] and studies the distribution of stress, strain, and wave propagation.However, these studies treat the powder as a continuum and determine the material constants experimentally.Force transmission in spherical particles occurs in a chain of contacts, which is usually referred to as the force chain.Force chains in granular matter have been widely investigated, experimentally [25,26] and in simulations [27,28].However these studies do not consider wave propagation velocity during loading.Liu and Nagel [29] and Jia et al. [30] found experimentally that sound propagates in granular media along strong force chains.Somfai et al. [31] investigate the sound waves propagation in a confined granular system.Recently Abd-Elhady et al. [32,33] studied contact time and force transferred due to an incident particle impact while using the Thornton and Ning approach [34] and found a good agreement between DEM simulation and experimental measurements.However, these studies are restricted to the particle collision and are not modeling the shock wave propagation during dynamic loading of particles.

Advances in Acoustics and Vibration
In the present work, DEM is used to simulate dynamic loading of a one-dimensional chain of spherical particles.The contact between particles is modeled using elastic and plastic loading, elastic unloading, and adhesion at contacts.Recently many researchers [35][36][37][38][39] used these contact models, however only to investigate different aspects of static compaction of particulate matter.In the current investigation, typical dynamic loading parameters are used, which are commonly found in high velocity compaction process.The 1D chain of spherical particles is chosen as a preliminary step towards the understanding of elastoplastic shock wave propagation and its effects during the entire loading process.Computer simulations reveal generation, transmission, and reflection of the elastoplastic shock wave through the particles.The shock wave effects on contact between particles, particle velocity, and its deformation are investigated.Effects of shock wave propagation on loading parameters are also investigated.In addition to transducer design, earthquake engineering, and soil mechanics, elastoplastic shock propagation in particulate materials has many other practical applications including the high-velocity compaction of powder material.

Basic Contact Equations
This section summarizes the theories that describe the contact behavior between particles and between particles and die wall during loading-unloading-reloading stages.Here the compact is modeled as a one-dimensional assembly of spherical particles that indent each other.These particles are assumed to be materially isotropic and homogeneous while depicting elastoplastic material behavior.Equivalent elastic modulus E * [38] is given by where E 0 represents Young's modulus and ν is Poisson's ratio of the material.The effective radius R 0 of the two particles in contact, here labeled 1 and 2 is determined by the relation As compaction proceeds, the particles overlap each other and elastic normal force follows the Hertzian law where h denotes the indentation or overlap between particles, as shown in Figure 1.In the plastic regime, as described by Storåkers et al. [40,41] for two spherical particles undergoing plastic deformation, the strain hardening relationship is given as where σ i is a material constant, M is the strain hardening exponent, and σ and ε are stress and strain in the uniaxial case.Normal contact force F p is given by the relation [42] F p = ηh (2+M)/2 , ( where Here σ 0 is a material parameter and, for ideally plastic material behavior, invariant c 2 = 1.43.By considering the material as perfectly plastic M = 0 while the particles having identical yielding stress normal contact force F p can be written as [39] The contact radius a is defined as [40] while contact stiffness The parameters F 0 , h 0 , and a 0 denote normal contact force, overlap, and contact radius, respectively, at the end of plastic compaction process, before the load is removed.The limit between a fully elastic unloading and a partially plastic unloading, as shown by Mesarovic and Johnson [42], is where p 0 is the uniform pressure at contact p 0 = 3σ y and w is the work of adhesion.In the present case, it is set to w = 0.5 J/m 2 .
During the unloading stage, the distance I 12 between the centers of the two particles which are pressed together is and the overlap where indentation recovered h u is given by [38,43] During unloading, the contact radius a is determined from ( 14) and normal force F u is given by [38,42] as which can be written in terms of χ as [42] In ( 15) and ( 16), the second term on the right hand side gives the contact force due to adhesion traction.During the present work it has been realized that the second term must be included in elastic and plastic loading equations when adhesion traction is considered.Both cases with and without adhesion traction are treated in the present study.

Discrete Element Method
The system under study is the one-dimensional dynamic compaction model shown in Figure 2.There is a chain of micron-sized identical, spherical particles aligned in a container with one end open and the other blocked.At the open end, these particles are in contact with a compaction tool which has the same diameter as those of the particles.Friction between the particles and the container walls is not considered.To start the compaction process, hydraulic pressure is used to accelerate the hammer which strikes the compact at a certain impact velocity.The hammer along with the compaction tool form the dynamic load.Compaction energy is mainly determined by the impact velocity and The dynamic compaction process is simulated by the discrete element method.This numerical method is used by Martin and Bouvard [37] and other authors to simulate static compaction.In the present work, DEM is used to extend fully developed contact models to simulate shock wave propagation in a chain of spherical particles.Here, each particle is modeled independently and interaction between neighboring particles is governed by contact laws as described in the previous section.This contact response plays an important role in the use of DEM to simulate shock wave propagation through the particles.During calculations at time t + Δt, where t is previous time and Δt is the time step, contact force between particles is calculated which determines the net force or compaction force F acting on each particle.By using Newton's second law, these resultant forces enable new acceleration, velocity, and position of each particle.At time t = 0, force, velocity, and position of each particle are known because it is the moment of the first hit.
The velocity v (t+Δt) i of a particle i at a time t+Δt is determined by adopting a central difference scheme as where the position x i is given by During iterative calculations, the size of time step Δt plays an important role to ensure numerical stability.For problems of a similar nature, Cundall and Strack [14] have proposed a relationship to calculate the time step which is further developed by O'Sullivan and Bray [44] for the central difference time integration scheme as where correction factor f t = 0.01 for the present case, m is the mass of the lightest particle, and k is the approximate contact stiffness given by expression (10).This value of the time step is shown to be sufficient to ensure numerical stability during the calculations.
During the compaction process, particle contact goes through several loading, unloading, and reloading sequences.In the beginning of compaction, contact force is initially elastic for small values of contact radius a and it is given by (3).The contact force follows the same curve during unloading and reloading in the elastic regime.At larger contact radius, the contact becomes plastic and contact force follows (8).The term relating the adhesion 2πwE * a 3/2 is added in elastic and plastic equations when considering adhesion traction.If the contact is unloaded, normal force follows elastic unloading (15).When contact is reloaded during unloading then it follows the same equation up to the value of the contact radius on which it was unloaded.Beyond this point, plasticity is reactivated and (8) applies.

Results and Discussion
For simulation, one hundred aluminum particles of diameter 2R = 100 μm are used.Dynamic compression load is applied on the particles by supplying a hydraulic pressure of 13.5 MPa which gives the hammer an impact velocity of 10 m/s.This impact velocity along with the different choices of loading mass results in a compaction energy of 1 J/g to 6.5 J/g.These loading parameters correspond to a typical high-velocity compaction process.The time step used is Δt = 2.3 ns which is estimated as explained in the previous section.The material properties of the aluminum particles are density 2700 kg/m 3 , yielding stress 146 MPa, Young's modulus 70 GPa, and Poissons's ratio 0.30.The material properties of the loading elements are density 7800 kg/m 3 , Young's modulus 210 GPa, and Poissons's ratio 0.35.The loading elements are made of steel with a high yielding stress, therefore, during loading, these elements deform only elastically.In the simulation, particle numbering starts from the compaction end.As a convention, resultant force, velocity,and displacement are taken positive from compaction to dead end, that is, along the positive x-axis, otherwise negative.Dynamic effects during particle compaction like elastoplastic shock wave propagation, particle contact behavior, and particle velocity along with loading parameters are investigated in this section.

Elastoplastic Wave Propagation.
The dynamic load transferred in particles is described using elastoplastic shock wave propagation variables like shock wave front velocity.The shock wave front is interpreted as the maximum absolute compaction force at a particular time while shock wave velocity is defined as the velocity of the wave front.The movement of the shock wave from compaction to dead end and then back to the compaction end is described as one compaction cycle.As the hammer moves forward to compact the material, particles overlap each other and thus contact forces are developed as a result of material stiffness and damping characteristics.The difference between contact forces results in a net force on the particle.This net force increases and the particle starts to move approximately with piston velocity after which this net force decreases and eventually becomes zero.During this period, the shock is also transferred continuously to the next particles.In one compaction cycle, the shock travels from the first particle to the last particle and then it is reflected back from the dead end towards the compaction end.During the backward part of the cycle, net contact force on a particle becomes negative as shown in Figure 3(a).The shock wave velocity is approximately 750 m/s (with adhesion) for the first cycle and it decreases approximately 5% from one cycle to the next as in Figure 3(a).Period in which shock passes through individual particle also changes from one cycle to the next and it is about 1.5 μs for the first half cycle.Figure 3(b) shows the enlarged view of neighboring particles.It can be seen that wave front is approximately shape preserving as it propagates through the particles.However, wave amplitude decreases slightly while shock wave passes from one particle to the next.It is mainly due to energy loss in the plastic deformation.In this particular case of a single chain of particles, shock  wave velocity and wave front mainly depend upon material properties and they are only slightly affected by changing the loading mass or initial impact velocity.

Particles Behavior during
Compaction.This section describes particle contact behavior, particle velocity, and compaction during the dynamic loading process.These parameters are mainly influenced by the shock wave propagation.
Contact history for different particles is shown in Figure 4(a).Contact response between neighboring particles plays an important role in transfer of mechanical energy through particles.The compaction process is initially elastic which remains for a very short time.Then particles are in plastic deformation where the overlap between particles increases linearly with contact force.However, contact force decreases at few points which depicts elastic unloading during compaction.It mainly occurs when shock wave travels back from the dead end and passes through the particle.Contact force increases as shock wave passes through the particle during both parts of the cycle, which can be seen in Figure 4(b).When hammer stops and there is no applied hydraulic pressure, contact force decreases and eventually becomes zero.During this elastic unloading, the small amount of overlap is recovered.Oscillations after unloading are due to adhesion traction between particles.In case of no adhesion between particles, those oscillations disappear.During compaction, velocity of the particles does not remain uniform and constant as illustrated in Figure 5.As the shock wave passes during forward part of the cycle, it results in the particle motion.Velocity of individual particle is increased from zero to approximately hammer velocity during this period.All particles start to move by the end of forward cycle after which shock hits the dead end and is reflected back.Now, as disturbance passes through the individual particle, its velocity decreases and eventually becomes zero.Time for motion of a particle is determined by the period between shock wave passes through the particle during forward and backward parts of the cycle.The hammer compresses particles until its velocity becomes zero.Up till this point, both cases of adhesion and no adhesion depict almost similar behavior.After compaction, net force becomes zero and particles expand and push back the hammer due to the elastic energy.In case of no adhesion, as in Figure 5(a), particles are moving with different velocities which indicates that particles are separated.In adhesion case, as illustrated in Figure 5(b), particles oscillations can be seen which are caused by adhesion between particles.It indicates that particles are not fully separated.Furthermore, shock wave propagation also plays an important role in particles deformation as shown in Figure 6(a).Particles are deformed plastically as shock passes during both parts of the cycle despite particles movement.However, displacement covered by the particles depends upon their position from the hammer as in Figure 6(b).All particles are compacted approximately to the same amount at shock wave propagation.

Effects of Changing Loading Parameters.
Like particles contact force and velocity, hammer kinetic energy (KE) is influenced by shock wave propagation.Piston KE and collective KE of all the particles are shown in Figure 7(a).There are three shock cycles for this compaction period.Hammer KE has the same pattern during a particular cycle and it changes when shock hits back the hammer at the end of the cycle.Particles KE increase as shock moves forward from compaction to the dead end.It reaches maximum value which is about 7% of hammer KE when shock hits the dead end.On the return cycle, particles KE decrease and become zero when shock hits the compaction end.At this point, all energy is converted to plastic deformation and elastic potential energy.For various choices of hammer KE, particles KE have different values but they all have the same pattern as illustrated in Figure 7(b).In all the cases, maximum and minimum value occurs when shock hits the dead end and compaction end, respectively.
It is obvious that particles compaction energy is mainly directly proportional to hammer kinetic energy which depends on its mass and impact velocity.Particles compaction for the the same hammer KE, but with different mass and velocity combinations is shown in Figure 8(a).
Compaction is the same in all the cases.However, compaction time is increased by increasing the loading mass.After compaction stage, loading mass oscillates with the same amplitude due to elastic energy of particles.However, oscillations period time is increased with increased hammer mass.In current simulation, the presence of hydraulic pressure after the compaction stage also served to avoid instantaneous particle separation.Contact force between the particles 52 and 53 for different values of hydraulic pressure is illustrated in Figure 8(b).It can be seen that oscillations are decreased with the increase in hydraulic pressure.In case the hydraulic force becomes larger than the contact force during unloading, then particles are compacted again plastically.

Conclusion
In the present numerical simulation, a discrete element method is employed to investigate elastoplastic shock propagation in a one-dimensional assembly of spherical particles.Well-established quasistatic compaction models are extended to the dynamic high-velocity range.Propagation and reflection of the elastoplastic shock wave in particles is simulated by using appropriate contact laws.Simulation results show that shock wave velocity, and shape of the wave front changed slightly during propagation.The shock wave determines the contact behavior, velocity and deformation of the particles during dynamic compaction.After compaction stage, adhesion traction restricts instantaneous particle separation.Particle deformation during one cycle initially remained almost the same regardless of loading parameters values.However, particles compaction depends upon kinetic energy of dynamic load despite different choices of impact velocity and loading mass.Shock wave propagation also affects the variations in hammer kinetic energy and collective kinetic energy of all the particles.Although the extension of the developed model into two and three dimensions requires more computational time and resources, it is nevertheless straightforward.

12 hFigure 1 :
Figure 1: Schematic of two particles in quasistatic contact.h is the overlap and I 12 is the distance between the centers of the two particles.

Figure 3 :
Figure 3: Elastoplastic shock wave during dynamic compaction.(a) Propagation and reflection during various cycles.(b) Shock front is shape preserving while amplitude decreases during propagation.

Figure 4 :Figure 5 :
Figure 4: Contact force between particles during loading-unloading process.(a) During compaction, contact is mainly in the plastic range except for elastic unloading at few points.(b) Contact force increases until the hammer stops.Oscillations during unloading are due to the adhesion between particles.

Figure 6 :
Figure 6: (a) Deformation of all the particles approximately remains the same.(b) Displacement covered by the particles depends upon their position from compaction end.

Figure 7 :
Figure 7: Kinetic energy varies as shock wave propagates through the particles.(a) The hammer kinetic energy and collective kinetic energy of all the particles.(b) Variation in collective kinetic energy of particles with hammer velocity.

Figure 8 :
Figure 8: (a) Compaction of the particles with the same hammer kinetic energy.(b) Hydraulic pressure counters the elastic unloading energy and also prevents the particles from instantaneous separation.