An algorithm for the grain-level modelling of a dry sand particulate system

In general, the recursive formula x_n+1 = R(x₁, x₂, ..., x_n) is commonly used to generate random numbers, in which R is a recursive function. According to the initial value (x₁, x₂,..., x_n), the new random number x_n+1 is derived. Obviously, the characteristic of this algorithm is that the random number sequence {x_n+1} is determined by the initial value (x₁, x₂,..., x_n) and the recursive function R. Therefore, the random number sequence {x_n+1} cannot meet the requirements of randomness and independence. The other characteristic of this algorithm is that when the sample of random numbers is large enough, the same random number sequence {x_n+1} may be obtained. Therefore, {x_n+1} is called a pseudo-random number. When the sample of random numbers is small, the pseudo-random number can meet the randomness completely.

In this research, the mixed congruent algorithm is used to generate pseudo-random numbers. The general form of the mixed congruent algorithm can be written as follows [36, 37]: if x_i is assumed to be the initial value, the recursive function is x_n+1 = N · x_n + C (mod M), in which N, C are constants, and M is a modulus operation.

Dry sand grains are mainly composed of crystalline silicon dioxide, i.e. quartz [1]. According to age, transport mechanism and mineral composition, sand grains have different shapes and sizes. Mahaney [7] revealed that most sand grains are random and convex in shape. In general, the shape of most old sand grains is random with rounded edges, and younger sand grains produced artificially by crushing sandstone have sharp edges and corners. Therefore, it is reasonable and realistic to model sand grains using 3D random convex polyhedra.

In this section, we propose an algorithm for developing a 3D random polyhedron to model sand grains. The developed model is an irregular polyhedron, random in shape and size. The random polyhedron is a convex solid with an arbitrary number of boundary surfaces, which is varied from 12 to 18 in this paper. In fact, the number of surfaces of a real sand grain may be larger than 18 or less than 12. Statistically, it is acceptable to use surfaces numbering from 12 to 18 in simulation.

The generation algorithm for sand grains is composed of two steps. The first is to generate a 3D random octahedron. The second step is to generate a 3D random polyhedron. The first step is described in detail as follows.

Step 1. The algorithm to generate a random quadrilateral, Q₁

Sub-step 1. According to the size, S_sd, of sand grains, generate a circle with radius S_sd/2.

Sub-step 2. Randomly generate four points S₁, S₂, S₃ and S₄ on the circle.

Sub-step 3. Generate a quadrilateral connecting the four vertices in turn. L₁, L₂, L₃, L₄ and β₁, β₂, β₃, β₄ are the lengths of the four edges, and the internal angles of the quadrilateral, respectively (figure 1).

Zoom In Zoom Out Reset image size

The above algorithm, which is referred to as Q₁, is a very simple way to generate a random cyclic quadrilateral. However, we find that it is not robust, as illustrated in figure 2. The quadrilateral in figure 2 is very sharp with very small angle β_i and very short edge L_i(i = 1, 4). The relationship between internal angle β_i and side length L_i of the quadrilateral can be expressed as

$$\begin{equation} \left. {\begin{array}{@{}l@{}} \cos \beta_{1} =\frac{L_{1}^{2} +L_{4}^{2} -L_{2}^{2} -L_{3}^{2}}{2(L_{1} \cdot L_{4} +L_{2} \cdot L_{3} )} \\ \cos \beta_{2} =\frac{L_{1}^{2} +L_{2}^{2} -L_{3}^{2} -L_{4}^{2}}{2(L_{1} \cdot L_{2} +L_{3} \cdot L_{4} )} \\ \cos \beta_{3} =\frac{L_{2}^{2} +L_{3}^{2} -L_{4}^{2} -L_{1}^{2}}{2(L_{2} \cdot L_{3} +L_{4} \cdot L_{1} )} \\ \cos \beta_{4} =\frac{L_{4}^{2} +L_{1}^{2} -L_{3}^{2} -L_{2}^{2}}{2(L_{4} \cdot L_{1} +L_{3} \cdot L_{2} )} \\ \end{array}} \right\}. \end{equation} \tag{ 1 }$$

From the above equations it can be seen that β_i can be controlled by L_i(i = 1, 4). The side length, L_i, can be written as follows:

$$\begin{equation} L_{i} =\left| {{\begin{array}{*{20}c} {x_{i+1} -x_{i}} \\ {y_{i+1} -y_{i}} \\ {z_{i+1} -z_{i}} \\ \end{array}}} \right|,\quad (i = 1, 4), \end{equation} \tag{ 2 }$$

where x₁ = x₅, y₁ = y₅, z₁ = z₅.

Zoom In Zoom Out Reset image size

The following relationship should be incorporated into the algorithm Q₁ to control the random cyclic quadrilateral shape.

$$\begin{equation} L_{i} >S_{{\rm sd}} /\psi_{{\rm 1}} (i={\rm 1},{\rm 4}), \end{equation} \tag{ 3 }$$

in which ψ₁ is the controlling parameter and usually given a positive number. The effect of ψ₁ is to avoid the existence of small internal angles and short edges. ψ₁ is set to be 3.0 in this paper.

Step 2. The algorithm to generate a random convex polyhedron, Q2

First, generate random points, S₅ and S₆, as shown in figure 3. Then connect the points and vertices of the quadrilateral to generate a convex octahedral. The intersection point of the line S₅S₆ and the plane S₁S₂S₃S₄ is O. The length of |S₅O| and |S₆O| should satisfy the following relations:

$$\begin{equation} \left. {{\begin{array}{*{20}c} {S_{{\rm sd}} \cdot \psi_{2} >\left| {S_{5} O} \right|>S_{{\rm sd}} \cdot \psi_{3}} \\ {S_{{\rm sd}} \cdot \psi_{2} >\left| {S_{6} O} \right|>S_{{\rm sd}} \cdot \psi_{3}} \\ \end{array}}} \right\}. \end{equation} \tag{ 4 }$$

The controlling parameters, ψ₂ and ψ₃, are set to be 0.75 and 0.25, respectively, to control the shape of the octahedral.

Zoom In Zoom Out Reset image size

The second step is to make the octahedral grow into the convex polyhedron with a random shape. So, we proposed the vector dot-product method to judge the convex of the polyhedron. The detailed algorithm is shown as follows.

Sub-step 1. Find the longest edge S_iS_j of the generated octahedral. Plane i and plane j are the two surfaces attached to the edge of S_iS_j. $\overrightarrow {V_{i}}$ and $\overrightarrow V_{j}$ are the unit normal vectors to plane i and plane j, and a new vector can be obtained as follows:

$$\begin{equation} \left. {\begin{array}{l} \overrightarrow {V_{i}} =\left[ {{\begin{array}{*{20}c} {x_{i}} \\ {y_{i}} \\ {z_{i}} \\ \end{array}}} \right] \\ \overrightarrow {V_{j}} =\left[ {{\begin{array}{*{20}c} {x_{j}} \\ {y_{j}} \\ {z_{j}} \\ \end{array}}} \right] \\ \overrightarrow {V_{ij}} =\left[ {{\begin{array}{*{20}c} {x_{i} +x_{j}} \\ {y_{i} +y_{j}} \\ {z_{i} +z_{j}} \\ \end{array}}} \right] \\ \end{array}} \right\}, \end{equation} \tag{ 5 }$$

or in more simple form:

$$\begin{equation} \overrightarrow V_{ij} =\overrightarrow V_{i} +\overrightarrow V_{j}. \end{equation} \tag{ 6 }$$

Sub-step 2. Generate a random point S_k in the edge S_iS_j. |S_kS_j| and |S_kS_i| are the length of point S_k to point S_j and point S_i. The point S_k is the one that meets the requirements of not being close to the points S_i and S_j when the following relations are satisfied:

$$\begin{equation} \left. {{\begin{array}{*{20}c} {S_{{\rm sd}} \cdot \psi_{4} >\left| {S_{k} S_{i}} \right|>S_{{\rm sd}} \cdot \psi_{5}} \\ {S_{{\rm sd}} \cdot \psi_{5} >\left| {S_{k} S_{j}} \right|>S_{{\rm sd}} \cdot \psi_{6}} \\ {S_{{\rm sd}} \cdot \psi_{6} >\left| {S_{k} S_{n}} \right|>S_{{\rm sd}} \cdot \psi_{7}} \\ \end{array}}} \right\}, \end{equation} \tag{ 7 }$$

in which

$$\begin{equation} \left. {\begin{array}{l} x_{k} =x_{i} +(x_{j} -x_{i} )\cdot RN_{1\times 1} \\ y_{k} =y_{i} +(y_{j} -y_{i} )\cdot RN_{1\times 1} \\ z_{k} =z_{i} +(z_{j} -z_{i} )\cdot RN_{1\times 1} \\ \end{array}} \right\}, \end{equation} \tag{ 8 }$$

where RN_1×1 is a random number matrix generated in section 2.1.

$$\begin{equation} \left. {\begin{array}{l} \left| {S_{k} S_{j}} \right|=\left| {{\begin{array}{*{20}c} {x_{k} -x_{j}} \\ {y_{k} -y_{j}} \\ {z_{k} -z_{j}} \\ \end{array}}} \right| \\ \left| {S_{k} S_{i}} \right|=\left| {{\begin{array}{*{20}c} {x_{k} -x_{i}} \\ {y_{k} -y_{i}} \\ {z_{k} -z_{i}} \\ \end{array}}} \right| \\ \end{array}} \right\}, \end{equation} \tag{ 9 }$$

in which S_n is the newly generated random point from S_k along the direction $\overrightarrow {V_{ij}}$. The controlling parameters, ψ₄, ψ₅, ψ₆ and ψ₇, are set to be 0.4, 0.1, 0.25 and 0.05, respectively, to avoid sharp corners.

Sub-step 3. Generate a new random polyhedron, as shown in figure 4. Connect point S_nS_i, S_nS_j, S_nS_a, S_nS_b in turn to generate a new polyhedron.

Sub-step 4. Judge the convexity of the polyhedron.

Zoom In Zoom Out Reset image size

The vector dot-product method is introduced here to check the convexity of the polyhedron, as shown in figure 5. $\overrightarrow {V_{{\rm r}}}$ is the unit normal vector to surface S_aS_nS_j, given as

$$\begin{equation} \overrightarrow {V_{{\rm r}}} =\left[ {{\begin{array}{*{20}c} {x_{r}} \\ {y_{r}} \\ {z_{r}} \\ \end{array}}} \right]. \end{equation} \tag{ 10 }$$

$\overrightarrow {V_{{\rm ai}}}$ is the unit vector from point S_a to S_i, given as

$$\begin{equation} \overrightarrow {V_{{\rm ai}}} =\frac{1}{L_{{\rm ia}}}\left[ {{\begin{array}{*{20}c} {x_{{\rm i}} -x_{{\rm a}}} \\ {y_{{\rm i}} -y_{{\rm a}}} \\ {z_{{\rm i}} -z_{{\rm a}}} \\ \end{array}}} \right], \end{equation} \tag{ 11 }$$

where L_ia is the length of the edge S_iS_a.

Zoom In Zoom Out Reset image size

If the generated polyhedron is convex, all vertices, S_i(i = 1, m), should lie in the opposite side of the plane r, which can be described as follows:

$\overrightarrow {V_{r}} \cdot \overrightarrow {V_{{\rm ai}}} <0$, convex,

$\overrightarrow {V_{r}} \cdot \overrightarrow {V_{{\rm ai}}} >0$, non-convex.

In the following, the vector dot-product method is also used to check overlapping in section 2.3.

Sub-step 5. Repeat the steps 2, 3 and 4 until the number of boundary faces of the polyhedron reach the given number. In this paper, the number is in the range of 12 to 18. If the convexity conditions are not satisfied, the polyhedrons are discarded.

The typical generated 3 D models of sand grains are shown in figure 6.

Zoom In Zoom Out Reset image size

In this section, we employ the take-and-place algorithm, which is referred to as Q₃, to generate the initial dry sand model. This algorithm is composed of two sub-steps. The first step is to generate all the required sand grains according to the algorithm Q₁ and Q₂ in section 2.2. The second step is to first generate a random list for all generated sand grains and then take and place the generated sand grains one by one into the analytical model according to the above random list, and finally check the boundary conditions and overlap between the sand grains. Details of the algorithm Q₃ are described as follows.

Sub-step 1. Generate sand grains according to the algorithm Q₁ and Q₂ in section 2.2, in which the number of sand grains should be large enough to meet the following take-and-place requirement.

Sub-step 2. Generate a random list for each sand grain, preparing for the take-and- place stage. So, a set of random numbers, RN_n×3, is generated first.

$$\begin{equation} RN_{n\times {\rm 3}} =\left[ {{\begin{array}{*{20}c} {r_{11}} & {r_{12}} & {r_{13}} \\ \vdots & \vdots & \vdots \\ {r_{n1}} & {r_{n2}} & {r_{n3}} \\ \end{array}}} \right], \end{equation} \tag{ 12 }$$

in which 1 > r_ij > 0, i = 1, n, j = 1, 3.

Then calculate the random points in the domain of the specimen, RP.

$$\begin{equation*} RP=RN_{n\times 3} \times ([B]_{n} -[B]_{m} )+I_{n\times 3} \times [B]_{m}, \end{equation*}$$

in which

$$\begin{equation} \left. {\begin{array}{l} \left[ B \right]_{m} =\left[ {{\begin{array}{*{20}c} {x_{m}} \\ {y_{m}} \\ {z_{m}} \\ \end{array}}} \right] \\ \left[ B \right]_{n} =\left[ {{\begin{array}{*{20}c} {x_{n}} \\ {y_{n}} \\ {z_{n}} \\ \end{array}}} \right] \\ \end{array}} \right\}, \end{equation} \tag{ 13 }$$

where [B]_m and [B]_n are the model boundary coordinates, and x_n > x_m, y_n > y_m, z_n > z_m. I_n×3 is an n × 3 unit matrix, shown as follows:

$$\begin{equation} I_{n\times 3} =\left[ {{\begin{array}{*{20}c} 1 & 1 & 1 \\ \vdots & \vdots & \vdots \\ 1 & 1 & 1 \\ \end{array}}} \right]. \end{equation} \tag{ 14 }$$

Sub-step 3. Take and place sand grains one by one into the analytical model according to the random list for each sand grain RP.

In order to do so, first move the center points, O', to RP, and then check to ensure that all vertices satisfy the boundary conditions.

$$\begin{equation} SS_{n} =RP+SS_{o} \end{equation} \tag{ 15 }$$

$$\begin{equation} SS_{n} =\left[ {{\begin{array}{*{20}c} {S_{1}} & \cdots & {S_{n}} \\ \end{array}}} \right]^{{\rm T}}, \end{equation} \tag{ 16 }$$

$$\begin{equation} S_{i} =[{\begin{array}{*{20}c} {x_{i}} & {y_{i}} & {z_{i}} \\ \end{array}}] \end{equation} \tag{ 17 }$$

where SS_o and SS_n are vertices of grains before and after placing.

When the boundary conditions are not satisfied, take and place the next sand grain until the boundary conditions are satisfied. The boundary conditions are as follows:

$$\begin{equation} \left. {{\begin{array}{*{20}c} {x_{n} >x_{i} >x_{m}} \\ {y_{n} >y_{i} >y_{m}} \\ {z_{n} >z_{i} >z_{m}} \\ \end{array}}} \right\}, \end{equation} \tag{ 18 }$$

or simply:

$$\begin{equation*} [B]_{n} >S_{i} >[B]_{m}. \end{equation*}$$

Sub-step 4. Check the overlap with the grains that have already been placed in the model. If not satisfied, take and place the next grain until meeting the overlapping requirements.

The algorithm of the vector dot-product method introduced in section 2.2 is employed to check the overlap between the sand grains. When overlap exists between the grains I and J, the point k can be found which is located inside the grain J, and lies on the faces of the grain I. The vector P_m is the normal of the face m of the grain J. The vector P_k is vertical to the face m from the point k. The dot-product of P_m and P_k must be positive.

Sub-step 5. Calculate the volume of the sand grains in the model and the total number of the sand grains that have been taken and placed. If the number of sand grains reaches the requirement, stop the take-and-place process and go to sub-step 6. If not, go to sub-step 3 and repeat.

Sub-step 6. Output the data for the sand grains which are randomly distributed in the analytical model.

The process of boundary condition checking in sub-step 3 is quite time consuming. Also, the time complexity of overlap checking in sub-step 4 is O(N²). Therefore, we propose the random rotational method and limited moving method to improve the boundary condition checking and overlap checking, respectively, in order to reduce computational time, which will be described in detail as follows.

The random rotational method, which is referred to as Q₄, is introduced first.

SS_n is assumed to be the vertices of the newly placed sand grain. Move the origin of the coordinate system to the center point of grains, O'. Make the grain rotate randomly around O' with angles of α, β and γ around the x-axis, y-axis and z-axis. The rotational vertices, SS_r, is described as

$$\begin{eqnarray} \fl SS_{{\rm r}} = SS_{n} \left[ {{\begin{array}{@{}lll@{}} {\cos \beta \cos \gamma} & {\cos \beta \sin \gamma} & {{\rm -}\sin \beta} \\ {\sin \alpha \sin \beta \cos \gamma -\cos \alpha \sin \gamma} & {\sin \alpha \sin \beta \sin \gamma +\cos \alpha \cos \gamma} & {\sin \alpha \cos \beta} \\ {\cos \alpha \sin \beta \cos \gamma +\sin \alpha \sin \gamma} & {\cos \alpha \sin \beta \sin \gamma -\sin \alpha \cos \gamma} & {\cos \alpha \cos \beta} \\ \end{array}}} \right].\nonumber\\ \end{eqnarray} \tag{ 19 }$$

The limited moving method, which is referred to as Q₅, is described as follows.

As mentioned above, SS_n is the vertices of the newly placed sand grain. SS_m is the coordinates after moving. The relationship between SS_n and SS_m is

$$\begin{equation} SS_{m} =\psi_{7} RN_{n\times 3} +SS_{n}, \end{equation} \tag{ 20 }$$

in which ψ₇ is a controlling parameter and RN_n×3 is a matrix of a random number list.

The limited moving distance is described as

$$\begin{equation} (r_{i1}^{2} +r_{i2}^{2} +r_{i3}^{2} )\cdot \psi_{7}^{2} <L_{n} \quad (i = 1, n), \end{equation} \tag{ 21 }$$

where L_n is a positive number given artificially. In this paper, we set L_n as 0.1 * S_sd. If the moving distance is very large, the grain may overlap with the other particles. We think that the moving distance should be very little and L_n is set to be 0.1 * S_sd.

Rotate and move the newly placed sand grain in a limited range, simply when there is overlap with other grains. This can save a great deal of computational time. We optimize the placing algorithm in Q₃ by adding random rotation, and limited moving in boundary and overlap checking. Figure 7 shows the random distribution in space of sand grains which gives the initial uncompact model of the sand grains.

Zoom In Zoom Out Reset image size

In the process of 'take-and-place', the grains that do not satisfy the boundary conditions for more than 1000 times will be discarded.

Dry sand is a kind of porous matter, compactable in response to pressure. Porosity is an important indicator for dry sand. In this paper, porosity is defined as the voids volume percentage over the specimen total volume which is not filled by grains. Porosity is dependent on the size, shape and arrangement of sand grains. Also, porosity differs between consolidated and unconsolidated dry sand. Compaction may lead to rearrangement of sand grains, reducing the pores between sand grains and increasing contact area, bulk density and friction force. Therefore, it is necessary to develop the compacting step.

The DEM has been widely used in recent years, especially in the domain of granular matter research. In this section, we present a compaction algorithm to increase the bulk density of the initial model of the sand specimen generated in section 2.3, employing the DEM.

2.4.1. Motion equations.

In DEM, the motion of a particle is determined by Newton's second law. All particles are assumed to be rigid bodies. The translation equation for the i-th particle is as follows:

$$\begin{equation} m_{i} {\ddot{u}}_{i}(t)+C_{mi}\dot{u}_{i}(t)=F_{i}(u,t). \end{equation} \tag{ 22 }$$

The rotation equation for the i-th particle is as follows:

$$\begin{equation} \Gamma_{i} \ddot{\psi}_i(t)+C_{\Gamma i}\dot{\psi}_i(t)=M_i(u,t), \end{equation} \tag{ 23 }$$

where, F_i(u, t) and M_i(u, t) are the sum of all forces and torques applied on particle i; u_i and ψ_i denote the translational displacement and rotational angle of particle i; m_i and Γ_i stand for the mass matrix and rotational inertia matrix; C_mi and C_Γi are translational and rotational damping coefficients matrices.

In our study, the sum of all sand particles is N. While dropping, only gravity and contact forces act on a particle; thus F_i and M_i can be decomposed as

$$\begin{equation} F_{i} =M_{i} g+\mathop \sum\limits_{j=0,j\ne i}^{N-1} F_{ij} \end{equation} \tag{ 24 }$$

$$\begin{equation} M_{i} =\mathop \sum\limits_{j=0,j\ne i}^{N-1} R_{j} \cdot F_{ij} \end{equation} \tag{ 25 }$$

where M_ig and F_ij stand for the gravity and contact forces of particle i; R_j is a vector pointing from the center of mass of particle i to the contact point with particle j.

2.4.2. Contact model in DEM.

We assumed that the contact particles interact according to Coulomb's contact law with dry friction. The contact force, F_ij, between particle i and particle j, can be described as follows:

$$\begin{equation} F_{ij} =F_{ij}^{n} +F_{ij}^{s}, \end{equation} \tag{ 26 }$$

where $F_{ij}^{n}$ and are normal contact force and tangential friction force.

According to Coulomb's contact law, tangential friction force can be written as

$$\begin{equation} F_{ij}^{s} =\varepsilon \cdot \mu \cdot F_{ij}^{n} ~(0\le \varepsilon <1), \quad \hbox{for static} \end{equation} \tag{ 27 }$$

$$\begin{equation} F_{ij}^{s} =\mu \cdot F_{ij}^{n} ,\quad \hbox{for dynamic,} \end{equation} \tag{ 28 }$$

where μ is the friction coefficient.

The normal contact force, $F_{ij}^{n}$, is decomposed as follows:

$$\begin{equation} F_{ij}^{n} =f_{ij}^{ne} +f_{ij}^{nd} =k_{n} p_{ij}^{n} +k_{t} \ddot{u}, \end{equation} \tag{ 29 }$$

where $f_{ij}^{ne}$ and $f_{ij}^{nd}$ stand for the elastic and viscous parts of normal contact force; k_n and k_t are normal contact stiffness and damping coefficient; $p_{ij}^{n}$ and $\ddot{u}$ denote normal penetration displacement and normal velocity.

2.4.3. Contact detection.

In this section, we present the finite element grid according to the model developed above.

There is no overlapping but only contact between sand grains. Therefore, the sand grains can be meshed one by one, alone. The unstructured grid algorithm is adopted in generating the grid of sand grains. In this section, 3D tetrahedron mesh is used, employing the 3D constrained Delaunay algorithm [33, 35]. The following gives the algorithm for the grid in detail.

Sub-step 1. Move the origin of the coordinate system to the center point O' of the i-th sand grain to form a local coordinate system X'O'Y', as shown in figure 8. Initialize a cube large enough to encase the sand grain. Initialize an array of linked lists, Gr−Link(), combining a static structure (an array) and a dynamic structure (linked lists), as shown in figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Sub-step 2. Calculate the values of the vertices of all sand grains in local coordinates in X'O'Y' and store the value in Gr−Link(). Generate original aided grid, Gr-Ori, as shown in figure 10, in the initial cube belonging to the i−th sand grains.

Zoom In Zoom Out Reset image size

Sub-step 3. Insert all vertices of the i−th sand grains into the aided grid Gr-Ori according to Delaunay principles [33, 35]. The original grid is remeshed to form a new grid Gr-New. Update the original grid Gr-Ori and the array Gr−Link() every time a vertex is inserted, until the last one. Store the data of grid Gr-New.

Sub-step 4. Optimize the grid Gr-New by inserting random points. Generate the optimized grid data Gr_Opt and store in the array Gr−Link(). Remove the grid outside the sand grain. Update the data in Gr−Link().

Sub-step 5. Move the origin of the coordinate system to the center point of the next sand grain and repeat sub-step 1–4 until the last one.

Sub-step 6. Update the data in Gr−Link() and output the final grid Gr_out of the sand specimen.

The grid algorithm is programmed in Fortran; the FE grid can be implemented into such finite element codes as ANSYS, ABAQUS, AUTODYN and LS-DYNA. The grid of the sand particulate system is shown in figure 11.

Zoom In Zoom Out Reset image size

We employed the finite element code to model the nonlinear mechanical response of the sand particulate system. The hydrocode is a powerful and reliable approach and has been used extensively in modelling static and dynamic events. In this paper, the hydrocode LS-DYNA [39] is used. In a hydrocode, material models which properly describe the material behavior under loadings are necessary. The material model, HJC model [38], is used to model the nonlinear behavior of individual sand grains.

In the JHC material model, material type 111 in LS-DYNA, the pressure is expressed as a function of the volumetric strain including the effect of permanent crushing. The damage is accumulated as a function of the plastic volumetric strain, equivalent plastic strain and pressure. The normalized equivalent stress is defined as

$$\begin{equation} \sigma^{{\rm \ast}}=\frac{\sigma}{f_{{\rm c}}^{{\rm '}}}, \end{equation} \tag{ 30 }$$

where σ is the actual equivalent stress, and f_c' is the quasi-static uniaxial compressive strength. The expression is written as

$$\begin{equation} \sigma^{\ast}=[{\rm A}(1-D)+{\rm B}P^{\ast N}][1-{\rm C}\ln (\dot{{\varepsilon}}^{\ast})], \end{equation} \tag{ 31 }$$

where A, B and C are constants, D is the damage parameter, P^* = P/f_c' is the normalized pressure and $\dot{{\varepsilon}}^{\ast}=\dot{{\varepsilon}}/\dot{{\varepsilon}}_{0}$ is the dimensionless strain rate. The sand particle is modelled using the JHC material model.

The parameters for individual sand particles are as follows [38, 39]:

$$\begin{eqnarray*}&&\fl\rho ={\rm 265}0\,{\rm kg\,m}^{-3},\ P_{{\rm v}} ={\rm 31}.{\rm 3}\%,\ C_{{\rm p}} ={\rm 86}\,{\rm J\,kg}^{-1}\,{\rm K}, \ A =0.{\rm 79},\ B={\rm 1}.{\rm 6}, \ N =0.{\rm 61},\nonumber\\ &&\fl C =0.0{\rm 1},\ S_{{\rm max}} ={\rm 7}.0,\ D_{1} =0.0{\rm 4},\ D_{{\rm 2}} ={\rm 1}.0,\ P_{{\rm crush}} = 0.0{\rm 1},\ T= 0.0{\rm 44}\,{\rm GPa}. \end{eqnarray*}$$

The contact and friction effect between sand grains happens and consumes energy constantly under static and dynamic loadings. Therefore, this effect should be taken into account in the modelling. The automatic surface-to-surface contact [39] is used to model the contact and friction effects. In addition, the friction coefficient is defined as follows [39]:

$$\begin{equation} \mu_{{\rm c}} =\mu_{{\rm k}} +\left( {\mu_{{\rm s}} -\mu_{{\rm k}}} \right){\rm e}^{-{\rm DC}\cdot \left| {v_{{\rm rel}}} \right|}, \end{equation} \tag{ 32 }$$

where μ_s is the static friction coefficient, μ_k is the dynamic friction coefficient and v_rel is the relative velocity of the surfaces in contact. DC is exponential decay coefficient. In simulation, the static and dynamic friction coefficients among sand particles are set to be 0.3 according to the research [13, 40]. The 1-node integrated solid formulation of the tetrahedral element is adopted in the simulation of sand particles.

In this section, we employ the sand particulate system developed above and the FE code to study the uniaxial response of dry sand under static compression. In simulation, the sand specimen is compressed under uniaxial loading, shown in figure 12. The specimen is a cylinder, 5 cm in diameter and 8 cm in length. Lateral displacement of the specimen is constrained using a steel tube. The specimen is constrained on the bottom surface. The quasi-static displacement loading is applied on the top surface. The computational parameters, material models and contact parameters in the FE code, are shown above. In the model, the steel tube is assumed to be a rigid body.

Zoom In Zoom Out Reset image size

The results are shown in figure 13; the axial stress and axial strain relationship. We find that there are three different stages in the uniaxial compression test. The first stage is approximately elastic compression. In this stage, the sand grains all experience elastic deformation. The friction effects between sand grains are not significant. However, this stage is very short. The maximum axial stress under elastic compression is up to 12.9 MPa, which is similar to test results [43]. Then, in the second stage, the specimen comes under high level compressive force. Obviously, the sand particulate system is compacted and the voids in the specimen decrease with higher axial compressive stress. Some sand grains experience a yielding stress state. The friction effects between sand grains are significant. We find a slight tumbling and rolling of the sand grains during compression in this stage. The maximum axial stress at this stage is up to 27.2 MPa, which is slightly less than the test results [43]. With the increase in pressure, the third stage starts. In this stage, the porosity increases because of compression leading to more surface contact. Most of the sand grains experience large deformation with plastic strain. The effects of contact and friction between the sand grains are very significant. Axial stress increases significantly. However, the axial strain increases slightly. This stage is the 'lock-up' stage [3], hardening the response. Figure 13 also shows that most grains experience damage at the third stage. The damage patterns forming along vertical bands in the three stages originate from the lateral confinement under the uniaxial loading. It reveals that the patterns of deformation and damage of sand grains are sensitive to the loading mode.

Zoom In Zoom Out Reset image size

Figure 14 gives the relationship between axial stress and void ratio under high axial pressures. When the void ratio ranges from 46.4% to 13.1%, the axial stress increases significantly, to the lateral confinement of the steel tube. The results agree well with the test data, in compression.

Zoom In Zoom Out Reset image size

In the simulation, the damage, D, of an individual sand grain is defined as $D=\sum\frac{\Delta \varepsilon_{{\rm p}} +\Delta \mu_{{\rm p}}}{D_{1} (P^{\ast}+T^{\ast})^{D_{2}}}$, where Δε_p and Δμ_p are the equivalent plastic strain and volumetric strain, D₁ and D₂ are defined above, and P^* and T^* are the normalized pressure and the normalized maximum tensile hydrostatic pressure. The damage of individual sand particles at different stages under uniaxial compression is shown in figure 13. We find that most sand grains experience elastic deformation at stage I. However, a high stress level leads to more and more damage to most sand particles at stage III. The damage of sand particles becomes more and more serious with the increase in stress level. The voids of the specimen become less and less. This leads to more contact surfaces between sand grains and higher axial stress, as shown in figure 13. Therefore, contact and friction play more important roles in stage III than in stage I. Individual sand particles are crushed by the contact force. Thus, the grain-level response, such as the deformation, damage and stress–strain response of individual sand particles, plays an important role and should be taken into account at a high stress level under static compression.

The dynamic response of sand grains is different to that of continuum media. Much experimental research [12, 41, 42, 49] has been carried out into the impact properties of dry sand. However, it is not yet well understood for dry sand under intense dynamic loadings. It is necessary to investigate the mechanical response of the dry sand particulate system under high rate loading originating from impact, gas and weapon explosions, earthquakes, etc. The study of the mechanical properties of sand particles under impact is very challenging. In this section, we present numerical approaches to study the grain level response of dry sand under impact. The numerical model of Split Hopkinson Pressure Bar (SHPB) is shown in figure 15. The computational parameters for sand particles are shown above. The bars and tube are assumed to be elastic with shear modulus G = 77 GPa and Poisson's ratio μ = 0.3.

Zoom In Zoom Out Reset image size

Kabir et al [8] carried out a series of tests on dry sand in 2010. First, we simulated the tests on the dry sand particulate system at a strain rate of 470 and 900 s⁻¹, to verify the reliability of the numerical analysis approach presented above. In the specimen, there are three different contact surfaces: between sand grains, between bars and the sand specimen, and between the tube and the sand specimen. The contact algorithm of 'automatic surface-to-surface' is used for these three contact surfaces. The dimensions of the sand specimen are 20 mm in diameter, 10 mm in thickness, 0.8 mm in diameter. The sand specimen is confined by a steel tube.

Figures 16 and 17 show the comparison of the results and the test data at strain rate 470 and 900 s⁻¹. It is shown that the numerical results of stress–strain curves agree well with the test data, indicating that the numerical approaches can give reliable predictions of the sand specimen under impact. It also can be seen that the simulation is softer than the test. However, there is a significantly stiffer response at around 10% strain in figure 16 and 18% strain in figure 17. This is because the porosity decreases with the increase in the axial strain. The lateral confinement plays a significant role in the axial stress of the numerical specimen at high strain rate loadings.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the simulation, sand particles are modelled using the JHC material model. The computational parameter C in the material model is a strain rate coefficient. When C is larger than 0, the strain rate effects of individual sand particle are conducted. In order to investigate the effects of strain rate sensitivity, we carried out simulations using different strain rate coefficients C, 0.0 and 0.09. As shown in figure 18, we can see that strain rate sensitivity of individual sand particle has little effect on the stress–strain curve under impact loading.

Zoom In Zoom Out Reset image size

In order to study the dynamic response of the sand particulate system, six studies are conducted considering different rate loadings at 1, 10, 100, 200, 470 and 900 s⁻¹. The effects of different rate loadings on stress–strain curves are shown in figure 19. It can be seen that the dynamic stress of the specimen increases with the increase in strain rate at the same strain. However, the stress value increases a little when strain rates vary from 1 to 900 s⁻¹. Therefore, the effects of the loading rate on dynamic response can be ignored. Tests [42, 49] also lead to the same conclusion. So, we can conclude that the dry sand particulate system is not significantly sensitive to loading rate.

Zoom In Zoom Out Reset image size

As is known to all, sand is an internal frictional material. The effects of friction between particles are studied in the following. The cases considering different friction coefficients, varying from 0.1 to 0.5, are simulated. We firstly simulate the frictionless case, setting the friction coefficients to be 0. All the cases are conducted at strain rate 470 s⁻¹. Figure 20 gives the results. The simulations reveal that the frictionless case reduces the specimen stress up to 40% more than other cases. Removing friction reduces the resistance of a sand specimen to external pressure, leading to large axial strain at the same stress level. When friction coefficients vary from 0.1 to 0.5, the stress level increases slightly. Moreover, the difference between the cases with friction coefficients 0.3 and 0.5 is less. The friction and non-friction cases show that the friction between sand is an important factor. Friction has significant effects on dynamic stress.

Zoom In Zoom Out Reset image size