|Computer Modeling in |
Engineering & Sciences
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: email@example.com
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.  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.  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.  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)  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 . The collision detection to establish if two polyhedra are in contact or not is called the GJK algorithm . The actual Minkowski overlap of two contacting polyhedra is obtained by the Expending Polytope Algorithm (EPA) . 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 , is defined by
The boundary of is denoted as . It is easy to establish that A and B is in contact or not is equivalent to that the origin is inside or outside . Fig. 1 shows three contact states of a quadrilateral (A) and a triangle (B) together with the corresponding Minkowski differences in relation to the origin.
When A and B are in contact, or the origin is inside , the contact (vector) distance between A and B is defined as the minimum distance that can be applied to B such that A and B are just in touch:
Then the contact overlap, denoted as , is the magnitude of , and is termed as the Minkowski overlap in , while the unit vector of defines the contact normal , i.e.,
The point on that has the shortest distance to the origin is called the contact point on . The two corresponding points of on the boundaries of A and B are respectively and , and they are called the contact points on A and B, respectively. A common contact point for both A and B can be taken as
Here, , , , , and are called the Minkowski contact features in , and they are illustrated in Fig. 2. Different contact types that are classified in  include vertex-edge, edge-edge for both polygons and polyhedra, and vertex-face and edge-face, face-face for polyhedra.
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 , while computing the Minkowski contact features is achieved by EPA  in the current work. Both algorithms are very effective and their details are also described in . The two initial issues associated with EPA for some special contact cases mentioned in  are resolved in .
2.3 Contact Force Models
After the Minkowski contact features of two polyhedra, mainly including overlap , normal and contact point , are obtained by EPA, the normal contact force can be evaluated as
where the function fn defines the magnitude of the force. As suggested in , fn can have a general power-law form
where kn is the stiffness coefficient, and is the exponent. A linear spring contact model is recovered when n = 1, while a “Hertz-type” contact model is obtained when n = 3/2. Although the contact model (4) seems to be the same as other overlap based contact models commonly used in DEM, the crucial difference for non-spherical particles is that in the current model both the overlap and the contact normal, together with the contact point where is applied, are the Minkowski contact features defined in a particular manner as described in Section 2.1. As formally proved in , using the Minkowski contact features in (4) for convex particles ensures that elastic energy will be conserved in any contact scenarios for an elastic impact. This is another energy-conserving normal contact model developed within the framework of the general contact theory established in , and is an alternative to the contact volume based contact models developed in  based on the same contact theory.
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; is the resultant contact force from all the contacts concerned; is the summation of all other applied forces on the particle; and and are the velocity and acceleration of the particle. The equation can be numerically solved used the standard central difference or leapfrog time integration scheme. No detail will be offered here.
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 and are the moment inertia tensor and angular velocity vector of the particle respectively, and represents the resultant moment/torque acting on the particle. An explicit damping term is omitted here, but can be included in . Note that , and are in the global coordinate system, and in particular is a function of time for non-spherical shapes. (7) can be expended into a form
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 diagonal matrix , and the corresponding three principal moment axes at the current time are denoted as a matrix in the global coordinate system. These three directions form a local coordinate system attached to the particle, with the mass centre of the particle as the origin. is related to by
As is a function of time when the particle rotates, so is . In the local coordinate system, (7) is reduced to the so-called Euler’s equation
where is the local angular velocity and is the local moment.
As is a diagonal matrix, (10) can be decoupled into component form
where is one of three permutations of . These decoupled equations can be numerically solved. Nevertheless, how to represent and update the orientation of the particle dictates the efficient and accuracy of solving the rotational motion of a non-spherical particle in DEM.
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 consists of one real component and three imaginary components and can be expressed by a four-element vector with one scalar plus a vector: . Geometrically a quaternion corresponds to a rotation state in space. In this work is assumed to be a unit quaternion:
A quaternion can be alternatively denoted in the following form which represents a rotation around a unit axis (vector) by an angle :
The inverse or conjugate of a quaternion is simply
which geometrically represents a rotation around the same axis with the same angle but in the opposite direction. Also .
Two rotations can be combined into one using quaternion multiplication. The multiplication of two quaternions and is defined as
which is equal to applying the rotation first and followed by the rotation . Quaternion multiplication is not commutative because the vector cross product is not commutative.
A quaternion can be applied to rotate a vector. The rotation of a vector by a quaternion is denoted as and defined by
where is the quaternion extension of the vector . The rotated vector assumes to take only the imaginary part of the quaternion .
Similarly, define an inverse rotation operator as
3.2 Numerical Time Integration of Angular Motion Using Quaternion
The time integration scheme presented below is based on , 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 with the time step . Assume that at time t = tn, and the orientation quaternion are obtained, and the total applied moment is given and remains constant in the interval. The task is to numerically solve Eq. (7) to obtain and at time t = tn+1.
First an intermediate angular moment is computed by
and then is transformed to the local coordinate system
In order to obtain the orientation , a series of incremental rotations around the three axes are required to be performed. Define an (incremental) quaternion associated with the (local) angular moment , the (local) moment inertia , a local rotation axis and the time step as
Further define an updated quaternion of a given quaternion in terms of , , and a direction i as
Now the orientation can be obtained from by applying five consecutive incremental rotations as
Finally the global angular velocity can be updated as
where is the local angular velocity which can be computed as
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 . 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:
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 is obtained by
where dmax is the maximum diameter, and is a random number between 0 and 1. The vectors from the vertices of each cell to the corresponding seed need to be shorten by this shrinking factor.
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.
4.2 Penetration Simulation
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 mm, mm and mm. The inner cavity of the projectile is not shown here but its effects on the mass and moments of inertia are properly considered. The projectile is also discretised into a polyhedron using the profile geometric parameters where 12 divisions are applied to the head part and 36 divisions are used in the circumference around the axis, so it can also be modelled by the same polyhedral DEM procedure. The composite structure target is not considered. The rockfill layer has the dimensions of .
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 10% of the mean size. The interaction between the projectile and the rockfill particles is visualised by the DEM simulation. The profile at the ZY plane of the penetration process for the 100 mm sample is shown in Fig. 7. It clearly shows that the affected range by the projectile is about 6–7 times of the average rock size with the target position at the centre.
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 to . As expected, the target with 290 mm rockfill has a very large deflection angle.
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.
|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.|