Computer Modeling in Engineering & Sciences |

DOI: 10.32604/cmes.2021.015913

ARTICLE

Discrete Element Modelling of Dynamic Behaviour of Rockfills for Resisting High Speed Projectile Penetration

1College of Mechanical and Vehical Engineering, Taiyuan University of Technology, Taiyuan, 030024, China

2Zienkiewicz Centre for Computational Engineering, Swansea University, Swansea, SA2 8PP, UK

*Corresponding Author: Y. T. Feng. Email: y.feng@swansea.ac.uk

Received: 26 January 2021; Accepted: 04 March 2021

Abstract: This paper presents a convex polyhedral based discrete element method for modelling the dynamic behaviour of rockfills for resisting high speed projectile penetration. The contact between two convex polyhedra is defined by the Minkowski overlap and determined by the GJK and EPA algorithm. The contact force is calculated by a Minkowski overlap based normal model. The rotational motion of polyhedral particles is solved by employing a quaternion based orientation representation scheme. The energy-conserving nature of the polyhedral DEM method ensures a robust and effective modelling of convex particle systems. The method is applied to simulate the dynamic behaviour of a rockfill system under impact of a high speed projectile. The rockfill sample is generated by a three-dimensional Voronoi meso method with a specific particle size distribution. The penetrating process of the projectile striking the rockfill target is simulated. Some physical quantities associated with the projectile such as the residual velocity, penetration resistance, and deflection angle are monitored which can reflect the influence of the characteristics of the rockfill target on its anti-penetration performance. It can be concluded that the developed polyhedral DEM method is a very promising numerical approach in analysing the dynamic behaviour of rockfill systems subject to high speed projectile impact.

Keywords: Discrete element method; minkowski overlap; polyhedral particles; rockfill protection system; high speed projectile penetration

In the protection engineering, rockfills have advantages of high strength, easy access and low cost which make them ideal protection materials to resist the penetration. A series of experimental studies on anti-penetration capabilities of various bursting layers consist of soil, rock or concrete, and bursting layers with different configurations are carried out by many scientific research groups [1–5]. The influence of configurations and contents of protection layers on the deflection of projectiles has been investigated which indicates that the yaw-inducing techniques are very effective for projectiles with a large ratio of length to diameter and the projectiles would be damaged by unsymmetrical impact force to some extents. Langheim et al. [2] carried out a series of penetration experiments with rockfill and concluded that good protection capability can be achieved by selecting proper size, gradation and compaction density of rockfill.

For the numerical simulation of rockfill protection problems, commercial software such as AUTODYN, LS-DYNA and ABAQUS are widely used. Mu et al. [6] conducted simulations of the rockfill target and concluded that a good protective capacity could be achieved when the size, thickness and strength of the rockfill target were matched with the projectile. Fang et al. [7] studied the penetration process of rockfill layers considering the random distribution of particle size and shape. In the study of penetration behavior of the reinforced concrete by Zhang et al. [8,9], the influence of size and fraction of coarse aggregates on the trajectory stability were considered using the meso-scale model.

As a rockfill system is inherently discrete and discontinuous, the Discrete Element Method (DEM) [10] can be used to simulate the interaction between rock particles and the projectile directly and obtain the energy dissipation, penetration depth and projectile velocity. Furthermore, the internal mechanism of macroscopic response can be revealed from the microscopic level. At present, the study of DEM modelling on dynamic behavior of the rock particle system under impact load is rare, and the traditional DEM with spherical elements is mostly adopted in some penetration problems [11,12]. The fundamental different behaviour of a spherical particle from a polyhedral particle makes the use of spheres to model rockfill systems less reliable.

The difficulty of using polyhedron as a basic discrete element lies in modelling of the contact behaviour between polyhedral particles. In this work, a newly developed Minkowski difference-based contact modelling methodology in [13,14] for convex polyhedral particles will be adopted to model the dynamic behaviour of rockfill systems under projectile impact. In this methodology, the contact overlap between two polyhedra is taken to be their Minkowski difference distance, or simply the Minkowski overlap [14]. The collision detection to establish if two polyhedra are in contact or not is called the GJK algorithm [15]. The actual Minkowski overlap of two contacting polyhedra is obtained by the Expending Polytope Algorithm (EPA) [16]. The contact force is computed based on a linear normal contact model. Using the Minkowski overlap and related contact features in the normal model ensures that the total elastic energy will be conserved for convex polyhedra in an elastic impact, which leads to robust discrete element simulations of complex problems including rockfill systems.

The paper is organised as follows. The next section will briefly introduce all the aspects of a polyhedral DEM model based on the methodology developed in [13,14], including the definition of Minkowski overlap and related contact features, the procedures to compute the features and the contact models used. Section 3 describes the quaternion based orientation representation of polyhedra and a symplectic time integration scheme to solve the angular motion of polyhedral particles. Numerical simulations of a rockfill system are conducted in Section 4 where a rock particle generation scheme is described first, followed by the description of all the model and parameters used for the simulations. Effects of the system parameters on the dynamic behaviour of the rockfill system subject to a high speed projectile impact are presented and discussed. Some conclusions are drawn in Section 5.

2 Convex Polyhedral Discrete Element Model

Consider two convex polyhedral particles in contact. Each polyhedron can have any number of vertices and flat faces. Each face is a convex polygon which can also have any number of vertices/edges. The discrete element modelling of the contact between these two polyhedra should establish if they are in contact, and if so then completely defines the contact features, including contact overlap (or volume), normal direction, contact point and contact force. This section will briefly describe the relevant definitions, computational procedures and contact models used for polyhedra in the current work, which are collectively called the polyhedral DEM model. The description below is mainly adopted from [13,14] where more details can be found.

Note that the polyhedral DEM model presented below is equally applicable to convex polygons. For simplicity, convex polygons will be used for all illustrations.

2.1 Minkowski Contact Features

Let A and B be the two polyhedra concerned. The Minkowski difference of A and B, denoted as

The boundary of

When A and B are in contact, or the origin is inside

Then the contact overlap, denoted as

The point

Here,

2.2 Contact Detection and Evaluation of Contact Features

The whole contact detection procedure in DEM includes three steps: Global search, local resolution and contact force calculation. In the first step, particles concerned, regardless of their actual shapes, are approximated by a simple enclosure (a bounding sphere or box). The objective of the global search is to establish a potential contact list for each particle. In the second step, the contact between two particles in the contact list is further checked based on their actual geometric shapes to establish the real contact state: separation, or overlap. For a pair of particles which are in contact, their contact forces will be evaluated in the third step.

As many effective global search algorithms already exist, there will be no further discussion about the global search. The third force evaluation step will be described in the next subsection. This subsection will focus on the local contact resolution for two polyhedra. The computational procedures required to fulfil the objective are to effectively determine if two polyhedra are in contact or not, and if so to obtain the Minkowski contact features defined in the preceding subsection.

Determining if two convex polyhedra are in contact or not is conducted by GJK [15], while computing the Minkowski contact features is achieved by EPA [16] in the current work. Both algorithms are very effective and their details are also described in [13]. The two initial issues associated with EPA for some special contact cases mentioned in [13] are resolved in [14].

After the Minkowski contact features of two polyhedra, mainly including overlap

where the function fn defines the magnitude of the force. As suggested in [14], fn can have a general power-law form

where kn is the stiffness coefficient, and

The tangential contact will be modelled by the classic Coulomb friction law. No detail needs to be described.

3 Rotational Motion of Polyhedral Particles

When all contact forces are evaluated for a particle, the dynamic governing equation for the translational motion of the particle can be expressed as

where m is the mass of the particle; c is the coefficient of viscous damping;

It is the angular motion of non-spherical particles that imposes significant challenges to discrete element modelling. Based on Newton’s second law, the governing equation for the rotational motion of a particle can be written as

where

This is clearly a highly non-linear second order differential equation, and imposes some challenges to obtain its numerical solution accurately and effectively.

For a given non-spherical particle in general, its three principal moments of inertia are assumed available and denoted by a

As

where

As

where

3.1 Orientation Representation of Polyhedral Particles

Unlike the centrally symmetric spherical particle, representing the orientation of a polyhedral particle is critical to the polyhedral DEM model for the determination of rotational motion and contact detection. In this work, the orientation of a polyhedron is represented by a quaternion which provides the most effective way to handle 3D rotation of a non-spherical shape.

A quaternion

A quaternion can be alternatively denoted in the following form which represents a rotation around a unit axis (vector)

The inverse or conjugate of a quaternion

which geometrically represents a rotation around the same axis

Two rotations can be combined into one using quaternion multiplication. The multiplication of two quaternions

which is equal to applying the rotation

A quaternion can be applied to rotate a vector. The rotation of a vector

where

Similarly, define an inverse rotation operator

3.2 Numerical Time Integration of Angular Motion Using Quaternion

The time integration scheme presented below is based on [19], which is one of symplectic integers that possess some superb long term conservation properties and therefore are widely used in general non-linear dynamics.

Consider a discretised time interval

First an intermediate angular moment

and then is transformed to the local coordinate system

In order to obtain the orientation

where

Further define an updated quaternion

Now the orientation

Finally the global angular velocity

where

4 Numerical Simulation of a Rockfill System

4.1 Generation of the Rockfill Sample

The polyhedral DEM outlined in the previous sections has the advantage of modelling irregular convex polyhedral particles directly. A rockfill system consisting of polyhedral rock particles is generated through the three-dimensional Voronoi meso model [8]. A graded rockfill sample can be obtained by introducing random parameters for shape and size control.

Firstly, N seeds are randomly placed into a given region V. Then the region is divided into N cells using the Voronoi tessellation technique, following the principle that the distance of arbitrary point within the polyhedron from its corresponding seed Si is less than the distance from any other seed. A vertex P of the polyhedron satisfies the following equation:

in which

The region volume V, the number of cells N and the average effective radius of cells r (i.e., the average distance between two seeds) have the following relation:

The initial 3D Voronoi diagram is shown in Fig. 3, in which each cell is a convex polyhedron that can be considered as a rockfill particle. The seamless cells need to be shrunk to obtain the scattered rocks. The shrinking factor of vertices within one cell must be the same to keep the convexity of the polyhedron. The shrinking factors for different cells are determined according to the actual grade of the rockfill sample. The shrinking factor for polyhedra with the size in the range

where dmax is the maximum diameter, and

The vertices and connectivities of the polyhedra obtained from the above procedure are imported to our polyhedral DEM procedure. However, the polyhedral sample generated by the Voronoi technique is very loose as shown in Fig. 4. To obtain a condense rockfill target, all rockfills drop under gravity and then are compressed by moving the top wall boundary downwards until a desired packing configuration is achieved.

The projectile penetrating the rockfill layer of a composite structure target is simulated using the above polyhedral DEM. The projectile is an axi-symmetric object with its cross-section profile shown in Fig. 5, with

The penetration model of the rockfill layer is generated by the method described in Section 4.1 and displayed in Fig. 6 with the projectile. The projectile strikes the rockfill layer from the centre of the left face at around 0.5 ms with an initial velocity of 1200 m/s and completely penetrates through the layer at around 3.5 ms.

A number of microscopic parameters need to be selected to describe the contact behaviour between the rock particles and also between the particles and the projectile. These include the normal stiffness kn of the linear normal contact model, the tangential stiffness ks, the coefficient of friction f and the coefficient of restitution e. The sensitive analysis of the parameters shows that the stiffnesses and the coefficient of friction have a minor influence on the penetration results, while the coefficient of restitution affects the final results significantly.

The velocity of the projectile decreases rapidly which cannot penetrate the rockfill layer when a small restitution is selected between the projectile and the particles. It is difficult to measure the trajectory and penetration resistance of the projectile in the field experiment. Therefore, the microscopic parameters cannot be determined precisely by adjusting the values to make the simulation results consistent with the experimental results. At the current stage, the main purpose of the penetration simulating is to show the potential benefit of modelling the interaction between the projectile and the rockfill using the polyhedral DEM. The values for all the parameters are selected as shown in Tab. 1, where er and eb are the coefficient of restitution between the rock particles and between the particles and the projectile respectively. With these parameters, the residual velocity of the projectile is almost the same as the experimental results.

A series of penetration simulations has been conducted using the above polyhedral DEM model and microscopic parameters. The mean size of rockfill particles is selected as 100, 130, 175, 200, 290 mm, respectively. The particle size obeys a normal distribution within the range of

The residual velocity, the penetration resistance and the deflection angle curves are shown in Fig. 8. The five residual velocity curves locate in two regions. For the four curves of the rockfill size between 100 to 200 mm, the residual velocity is from 1000 to 1025 m/s. The residual velocity of the 290 mm rockfill target is much larger than others as 1100 m/s. As the diameter of the projectile is around 150 mm, we can deduce that the target whose rockfill size is closer to the projectile diameter has a better velocity reduction effect. For the penetration resistance, the fluctuation range of the sample with the large size rockfill is larger than that of the sample with the small size rockfill (100 and 130 mm). It is because the resistance is the resultant force of all the rockfill acting on the projectile, the number of rockfill contacted with the projectile is more stable for the target with a small size rockfill. The curves of the deflection angle show a similar tendency with the curves of the residual velocity. The deflection angle for all the targets stays in the range from

This paper has presented a polyhedral DEM method for modelling the dynamic behaviour of rockfills for resisting high speed projectile penetration. The key issues of the polyhedral DEM method have been systematically introduced including the calculation of the contact force and rotational motion of polyhedral particles. To determine the contact force, the contact relation of two convex polyhedra is conducted by GJK and then the Minkowski contact features are computed by EPA. This normal contact model is energy-conserving for any contact scenario of an elastic impact. To calculate the rotational motion, the orientation of polyhedral particles is expressed by a quaternion and the time integration is based on a symplectic splitting method. This polyhedral DEM method has been applied to model the dynamic behaviour of a rockfill system. The rockfill sample is generated by a three-dimensional Voronoi meso model with a specific particle size distribution. The penetrating process of a high speed projectile striking a rockfill target has been simulated. Some process variables of the projectile such as the residual velocity, penetration resistance, and deflection angle is monitored which can reflect the influence of the characteristics of the rockfill target on its anti-penetration performance. It can be concluded that the developed polyhedral DEM method is a potentially powerful numerical approach in analysing the dynamic behaviour of rockfill systems.

Acknowledgement: The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.

Funding Statement: This work is partially supported by National Natural Science Foundation of China under Grant No. 12072217 and by Open Fund of State Key Laboratory of Coal Resources and Safe Mining, China University of Mining and Technology, Beijing, China [Grant No. SKLCRSM19KFA12]. The support is gratefully acknowledged.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

- Austim, C., Halsey, C., Clodt, R. (1982). Protection systems development. Florida, USA: Tyndall Air Force Base.
- Langheim, H., Schmolinske, E., Stilp, A., Pahl, H. (1993). Subscale penetration tests with bombs and advanced penetrators against hardened structures. 6th International Symposium, Interaction of Nonnuclear Munitions with Structures, Florida, USA.
- Gelman, M., Richard, B., Ito, Y. (1987). Impact of armor-piercing projectile into array of large caliber boulders. Mississippi, USA: Waterways Experiment Station.
- Rohani, B. (1987). Penetration of kinetic energy projectiles into rock-bubble/boulder overlays. Mississippi, USA: Waterways Experiment Station.
- Underwood, J. (1995). Effectiveness of yaw-inducing deflection grids for defeating advanced penetrating weapons. Technical report. Panama City FL: Applied Research Associates, Inc.
- Mu, C., Shi, P., & Xin, K. (2012). Numerical simulation on rock anti-penetration layer against penetrating.
*Ordnance Material Science and Engineering*,*35(5)*, 4-8. [Google Scholar] - Fang, Q., & Zhang, J. (2014). 3D numerical modeling of projectile penetration into rock-rubble overlays accounting for random distribution of rock-rubble.
*International Journal of Impact Engineering*,*63(5)*, 118-128. [Google Scholar] [CrossRef] - Zhang, J., Yang, H., Zhang, Y., Wang, Z., & Wang, Z. (2018). Numerical investigations on the effect of reinforcement on penetration resistance of concrete slabs using a 3D meso-scale method.
*Construction and Building Materials*,*188*, 793-808. [Google Scholar] [CrossRef] - Zhang, J., Chen, W., Hao, H., Wang, Z., & Wang, Z. (2020). Performance of concrete targets mixed with coarse aggregates against rigid projectile impact.
*International Journal of Impact Engineering*,*141(2)*, 103565. [Google Scholar] [CrossRef] - Cundall, P. A., & Strack, O. D. (1979). A discrete numerical model for granular assemblies.
*Geotechnique*,*29(1)*, 47-65. [Google Scholar] [CrossRef] - Kudryavtsev, O., & Sapozhnikov, S. (2016). Numerical simulations of ceramic target subjected to ballistic impact using combined DEM/FEM approach.
*International Journal of Mechanical Sciences*,*114(3)*, 60-70. [Google Scholar] [CrossRef] - Holmen, J. K., Olovsson, L., & BĂ¸rvik, T. (2017). Discrete modeling of low-velocity penetration in sand.
*Computers and Geotechnics*,*86*, 21-32. [Google Scholar] [CrossRef] - Feng, Y., & Tan, Y. (2020). On minkowski difference-based contact detection in discrete/discontinuous modelling of convex polygons/polyhedra.
*Engineering Computations*,*37(1)*, 54-72. [Google Scholar] - Feng, Y., Tan, Y. (2021). The minkowski overlap and the energy conserving contact model for discrete element modelling of convex non-spherical particles. International Journal for Numerical Methods in Engineering.
- Gilbert, E. G., Johnson, D. W., & Keerthi, S. S. (1988). A fast procedure for computing the distance between complex objects in three-dimensional space.
*IEEE Journal on Robotics and Automation*,*4(2)*, 193-203. [Google Scholar] [CrossRef] - van den Bergen, G. (2001). Proximity queries and penetration depth computation on 3D game objects.
*Game Developers Conference:**Vol. 170*,San Jose, USA. - Feng, Y. (2021). An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles: Basic framework and general contact model.
*Computer Methods in Applied Mechanics and Engineering*,*373(1)*, 113454. [Google Scholar] [CrossRef] - Feng, Y. (2021). An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles: Contact volume based model and computational issues.
*Computer Methods in Applied Mechanics and Engineering*,*373(3)*, 113493. [Google Scholar] [CrossRef] - Dullweber, A., Leimkuhler, B., & McLachlan, R. (1997). Symplectic splitting methods for rigid body molecular dynamics.
*Journal of Chemical Physics*,*107(15)*, 5840-5851. [Google Scholar] [CrossRef]

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. |