- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Article Processing Charges ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Recently Accepted Articles ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
Advances in Acoustics and Vibration
Volume 2011 (2011), Article ID 123695, 9 pages
Discrete Element Simulation of Elastoplastic Shock Wave Propagation in Spherical Particles
Marcus Wallenberg Laboratory for Sound and Vibration Research (MWL), Royal Institute of Technology (KTH), 100 44 Stockholm, Sweden
Received 18 October 2010; Accepted 15 June 2011
Academic Editor: Mohammad Tawfik
Copyright © 2011 M. Shoaib and L. Kari. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Elastoplastic shock wave propagation in a one-dimensional assembly of spherical metal particles is presented by extending well-established 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.
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–3]. The dynamic wave propagation in granular media shows distinct behavior from the wave propagation in continues media . Shukla and Damania  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  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.  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–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.  perform numerical simulations to investigate the effects of the contact laws on wave propagation in granular matter. Similarly Sadd et al.  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  and this numerical method has been widely used for granular material simulations [15–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–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  and Jia et al.  found experimentally that sound propagates in granular media along strong force chains. Somfai et al.  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  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.
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–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.
2. 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  is given by where represents Young’s modulus and is Poisson’s ratio of the material. The effective radius 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 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 is a material constant, is the strain hardening exponent, and and are stress and strain in the uniaxial case. Normal contact force is given by the relation  where Here is a material parameter and, for ideally plastic material behavior, invariant . By considering the material as perfectly plastic while the particles having identical yielding stress normal contact force can be written as  The contact radius is defined as  while contact stiffness
The parameters , , and 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 , is where is the uniform pressure at contact and is the work of adhesion. In the present case, it is set to J/m2.
During unloading, the contact radius is determined from (14) and normal force is given by [38, 42] as which can be written in terms of as  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.
3. 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 mass of the dynamic load. Hydraulic pressure is maintained during the complete compaction process.
The dynamic compaction process is simulated by the discrete element method. This numerical method is used by Martin and Bouvard  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 , where is previous time and is the time step, contact force between particles is calculated which determines the net force or compaction force acting on each particle. By using Newton’s second law, these resultant forces enable new acceleration, velocity, and position of each particle. At time , force, velocity, and position of each particle are known because it is the moment of the first hit. The velocity of a particle at a time is determined by adopting a central difference scheme as where the position is given by
During iterative calculations, the size of time step plays an important role to ensure numerical stability. For problems of a similar nature, Cundall and Strack  have proposed a relationship to calculate the time step which is further developed by O’Sullivan and Bray  for the central difference time integration scheme as where correction factor for the present case, is the mass of the lightest particle, and 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 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 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.
4. Results and Discussion
For simulation, one hundred aluminum particles of diameter = 100 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 ns which is estimated as explained in the previous section. The material properties of the aluminum particles are density 2700 kg/m3, 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/m3, 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 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.
4.1. 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.
4.2. 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.
4.3. 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.
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.
The authors gratefully acknowledge the Higher Education Commission (HEC) of Pakistan and the Swedish Institute (SI) for their financial contribution to the project. The Swedish Agency for Innovation Systems, VINNOVA, is also acknowleged for supporting this investigation under Contract no. 2004-02896 (VAMP 30).
- H. P. Rossmanith and A. Shukla, “Photoelastic investigation of dynamic load transfer in granular media,” Acta Mechanica, vol. 42, no. 3-4, pp. 211–225, 1982.
- C. H. Liu and S. R. Nagel, “Sound in a granular material: disorder and nonlinearity,” Physical Review B, vol. 48, no. 21, pp. 15646–15651, 1993.
- V. F. Nesterenko, Dynamics of Heterogeneous Materials, Springer, New York, NY, USA, 2001.
- A. Shukla and C. Damania, “Experimental investigation of wave velocity and dynamic contact stresses in an assembly of disks,” Experimental Mechanics, vol. 27, no. 3, pp. 268–281, 1987.
- A. Shukla and C. Y. Zhu, “Influence of the microstructure of granular media on wave propagation and dynamic load transfer,” Journal of Wave-Material Interaction, vol. 3, pp. 249–265, 1988.
- K. Tanaka, M. Nishida, T. Kunimochi, and T. Takagi, “Discrete element simulation and experiment for dynamic response of two-dimensional granular matter to the impact of a spherical projectile,” Powder Technology, vol. 124, no. 1-2, pp. 160–173, 2002.
- B. O. Hardin and G. E. Blandford, “Elasticity of particulate materials,” Journal of Geotechnical Engineering, vol. 115, no. 6, pp. 788–805, 1989.
- P. J. Digby, “Effective elastic moduli of porous granular rocks,” Journal of Applied Mechanics, vol. 48, no. 4, pp. 803–808, 1981.
- C. Thornton and D. J. Barnes, “Computer simulated deformation of compact granular assemblies,” Acta Mechanica, vol. 64, no. 1-2, pp. 45–61, 1986.
- C. S. Chang and L. Ma, “Modeling of discrete granulates as microploar continua,” Engineering Mechanics, vol. 116, pp. 2703–2721, 1990.
- K. Walton, “The effective elastic moduli of a random packing of spheres,” Journal of the Mechanics and Physics of Solids, vol. 35, no. 2, pp. 213–226, 1987.
- M. H. Sadd, Q. Tai, and A. Shukla, “Contact law effects on wave propagation in particulate materials using distinct element modeling,” International Journal of Non-Linear Mechanics, vol. 28, no. 2, pp. 251–265, 1993.
- M. H. Sadd, G. Adhikari, and F. Cardoso, “DEM simulation of wave propagation in granular materials,” Powder Technology, vol. 109, pp. 222–233, 2000.
- P. A. Cundall and O. D. L. Strack, “Discrete numerical model for granular assemblies,” Geotechnique, vol. 29, no. 1, pp. 47–65, 1979.
- Y. C. Zhou, B. D. Wright, R. Y. Yang, B. H. Xu, and A. B. Yu, “Rolling friction in the dynamic simulation of sandpile formation,” Physica A, vol. 269, no. 2, pp. 536–553, 1999.
- Y. F. Cheng, S. J. Guo, and H. Y. Lai, “Dynamic simulation of random packing of spherical particles,” Powder Technology, vol. 107, no. 1-2, pp. 123–130, 2000.
- Z. P. Zhang, L. F. Liu, Y. D. Yuan, and A. B. Yu, “A simulation study of the effects of dynamic variables on the packing of spheres,” Powder Technology, vol. 116, no. 1, pp. 23–32, 2001.
- I. Goldhirsch and C. Goldenberg, “Stress in dense granular materials,” in The Physics of Granular Media, pp. 3–22, Wiley-VCH, 2004.
- J. Duran, Sand Powders and Grains, An Introduction to the Physics of Granular Materials, Springer, Berlin, Germany, 2000.
- N. Sukegawa, T. Sano, S. Horikoshi, and H. Takeishi, “Dynamic powder compaction for parts with high-aspect ratio,” International Journal of Impact Engineering, vol. 24, no. 6, pp. 561–570, 2000.
- R. S. Ransing, D. T. Gethin, A. R. Khoei, P. Mosbah, and R. W. Lewis, “Powder compaction modelling via the discrete and finite element method,” Materials and Design, vol. 21, no. 4, pp. 263–269, 2000.
- A. F. Bower, N. A. Fleck, A. Needleman, and N. Ogbonna, “Indentation of a power law creeping solid,” in Proceedings of the Royal Society of London, vol. 441, pp. 97–124, 1993.
- B. Storåkers and P. L. Larsson, “On Brinell and Boussinesq indentation of creeping solids,” Journal of the Mechanics and Physics of Solids, vol. 42, no. 2, pp. 307–332, 1994.
- N. Ogbonna, N. A. Fleck, and A. C. F. Cocks, “Transient creep analysis of ball indentation,” International Journal of Mechanical Sciences, vol. 37, no. 11, pp. 1179–1202, 1995.
- J. Geng, G. Reydellet, E. Clément, and R. P. Behringer, “Green's function measurements of force transmission in 2D granular materials,” Physica D, vol. 182, no. 3-4, pp. 274–303, 2003.
- D. W. Howell and R. P. Behringer, “Fluctuations in granular media,” Chaos, vol. 9, no. 3, pp. 559–572, 1999.
- F. Radjai, S. Roux, and J. J. Moreau, “Contact forces in a granular packing,” Chaos, vol. 9, no. 3, pp. 544–554, 1999.
- C. Goldenberg and I. Goldhirsch, “Force chains, microelasticity, and macroelasticity,” Physical Review Letters, vol. 89, no. 8, pp. 84302–084302, 2002.
- C. H. Liu and S. R. Nagel, “Sound in sand,” Physical Review Letters, vol. 68, no. 15, pp. 2301–2304, 1992.
- X. Jia, C. Caroli, and B. Velicky, “Ultrasound propagation in externally stressed granular media,” Physical Review Letters, vol. 82, no. 9, pp. 1863–1866, 1999.
- E. Somfai, J. N. Roux, J. H. Snoeijer, M. Van Hecke, and W. Van Saarloos, “Elastic wave propagation in confined granular systems,” Physical Review E, vol. 72, pp. 21301–21320, 2005.
- M. S. Abd-Elhady, C. C. M. Rindt, and A. A. van Steenhoven, “Contact time of an incident particle hitting a 2D bed of particles,” Powder Technology, vol. 191, no. 3, pp. 315–326, 2009.
- M. S. Abd-Elhady, C. C. M. Rindt, and A. A. V. Steenhoven, “Force propagation speed in a bed of particles due to an incident particle impact,” Advanced Powder Technology, vol. 21, no. 2, pp. 150–164, 2010.
- C. Thornton and Z. Ning, “A theoretical model for the stick/bounce behaviour of adhesive, elastic- plastic spheres,” Powder Technology, vol. 99, no. 2, pp. 154–162, 1998.
- O. Skrinjar and P. L. Larsson, “Cold compaction of composite powders with size ratio,” Acta Materialia, vol. 52, no. 7, pp. 1871–1884, 2004.
- O. Skrinjar and P. L. Larsson, “On discrete element modelling of compaction of powders with size ratio,” Computational Materials Science, vol. 31, no. 1-2, pp. 131–146, 2004.
- C. L. Martin and D. Bouvard, “Study of the cold compaction of composite powders by the discrete element method,” Acta Materialia, vol. 51, no. 2, pp. 373–386, 2003.
- C. L. Martin, “Unloading of powder compacts and their resulting tensile strength,” Acta Materialia, vol. 51, no. 15, pp. 4589–4602, 2003.
- C. L. Martin, “Elasticity, fracture and yielding of cold compacted metal powders,” Journal of the Mechanics and Physics of Solids, vol. 52, no. 8, pp. 1691–1717, 2004.
- B. Storåkers, S. Biwa, and P. L. Larsson, “Similarity analysis of inelastic contact,” International Journal of Solids and Structures, vol. 34, no. 24, pp. 3061–3083, 1997.
- B. Storåkers, N. A. Fleck, and R. M. McMeeking, “The viscoplastic compaction of composite powders,” Journal of the Mechanics and Physics of Solids, vol. 47, no. 4, pp. 785–815, 1999.
- S. D. Mesarovic and K. L. Johnson, “Adhesive contact of elastic-plastic spheres,” Journal of the Mechanics and Physics of Solids, vol. 48, no. 10, pp. 2009–2033, 2000.
- K. L. Johnson, Contact Mechanics, Cambridge University press, Cambridge, UK, 1985.
- C. O'Sullivan and J. D. Bray, “Selecting a suitable time step for discrete element simulations that use the central difference time integration scheme,” Engineering Computations, vol. 21, no. 2-4, pp. 278–303, 2004.