Abstract
A modified smoothed particle hydrodynamics (SPH) method is applied to simulate the failure process of heterogeneous materials. An elastoplastic damage model based on an extension form of the unified twin shear strength (UTSS) criterion is adopted. Polycrystalline modeling is introduced to generate the artificial microstructure of specimen for the dynamic simulation of Brazilian splitting test and uniaxial compression test. The strain rate effect on the predicted dynamic tensile and compressive strength is discussed. The final failure patterns and the dynamic strength increments demonstrate good agreements with experimental results. It is illustrated that the polycrystalline modeling approach combined with the SPH method is promising to simulate more complex failure process of heterogeneous materials.
1. Introduction
The microstructure of heterogeneous solids usually consists of several mineral components and some preexisting defects, such as void and cracks. With the advance of experimental techniques, close observations on the specimens’ activities towards their microstructures during the failure processes become possible. For instance, Tapponnier and Brace [1] investigated stress-induced microcrack developments within different mineral components in the Westerly granite. They noted that new cracks appeared firstly at grain boundaries. With a subsequent load, the transgranular cracks are developed. Around the peak stress, biotite failure including the kinking and sliding can be observed. They also suggested that biotite grains might control the strength of the granite. Wong [2] further investigated the faulting mechanisms of different minerals in Westerly granite with the confining pressure and temperature. It was concluded that the failure mechanisms were related to both mineralogy and the mineral grain orientation. Besides, several numerical models have been put forward to study the failure processes of heterogeneous materials, such as the lattice-based model [3], the RFPA code [4], the local degradation model based on the FLAC code [5], and the synthetic rock mass model (SRM) based on the particle flow code [6]. Even with great improvements, most of these models were based on finite element or discrete element methods with some limitations. For example, the finite element models cannot work well for the discontinuous failure process featured by distortion, large deformation, fracture propagation, and fragmentation. On the contrary, discrete element methods are not well defined in treating the continuous deformation and most of them cannot properly model the aggregates or grains in their microstructures as well as consider their effects in rock failure simulations.
This paper first introduces the SPH method and its governing equations. The developed SPH code can deal with both discontinuities and material heterogeneity especially for the simulation under dynamic loading conditions. An elastoplastic damage model is adopted to model the strength behavior of brittle solids. Then, a microstructure model is applied to generate the artificial granite specimen. A Voronoi diagram is firstly employed to discretise the physical domain; the generated polygons are classified into three types of minerals with different homogeneous index values for feldspar, quartz, and biotite. Simulations of the Brazilian splitting test and the uniaxial compression test choose two different granite samples based on the generated artificial specimen. The strain rate effects on the predicted dynamic tensile and compressive strengths are compared with those from experimental tests.
2. SPH Method and Governing Equations
The first application of the SPH method is to solve astrophysical problems [7, 8]. By combining the strength and damage model, this method was extended to simulate the fracture process of brittle solids [9]. The formulation of the SPH method is based on the pioneer work by Monaghan [10] and Libersky et al. [11]. Since the SPH method was more convenient than FEM on the simulation of fracture process, it was utilized to study the brittle failure and subsequent fragmentation in the damaged solids [12].
The foundation of SPH method is based on the interpolation theory. Assuming that is a kernel function which has a compact supporting domain with a radius of kh that gives the “kernel estimate” of the field variables at a point, the properties of each particle are integrated as sum values of its neighboring particles. An approximation form of a function and its differential form at point can be constructed as where the summation is over all the particles within the supporting domain of the given particle . These influenced particles are neighboring particles of the particle which have mass , position , density , velocity , and other properties. is the smoothing length and is the smoothing kernel function. The interpolation kernel called B-spline which has a compact support of 2 h is widely used in the SPH by Monaghan and Lattanzio [13], Libersky et al. [11], Randles and Libersky [14], and so on. This kernel interpolates to the second order in the smoothing length h and it is always nonnegative. The mass and momentum conservation equations of continuum mechanics are approximated as In the above equations, dependent variables include the density , the velocity , and the stress tensor . Independent variables contain spatial coordinates and the time . The total time derivative operator is acquired from the moving Lagrangian frame. The Greek superscripts and represent coordinate directions. The summation is taken over repeated Greek indices.
is calculated by an explicit time integration approach using the constitutive model through the strain-rate tensor which is expressed by the derivatives of the velocity as is the artificial viscous pressure aimed at smoothing shocks over a few resolution lengths and at stabilizing numerical solutions which are introduced by Gingold and Monaghan [8].
In this simulation, a neighbor search routine is implemented to list particles within their supporting domains. Once the interaction terms are computed, the evolution for each particle can be stepped forward in time by using the Leap-Frog integral method [8]. The Leap-Frog scheme is conditionally stable. As other explicit integral methods, its stability and accuracy can be guaranteed under the Courant-Friedrichs-Lewy (CFL) type condition. The time step is adopted by Libersky et al. [11] as where is the adiabatic sound speed for particle , is the particle speed, and is a constant factor around 0.3. In order to mitigate the tensile instability, the conservative smoothing method will be also performed to smooth the density, velocity, and energy fields [14].
3. Microstructure Modeling Method to Generate the Artificial Specimen
3.1. 2D Discretization of Domain Based on Voronoi Diagram
The discretization process takes a square domain as an example. The Voronoi diagram composed of many convex polygons is employed to explicitly discretise the physical domain. Those polygons are constructed by a set of spatial points filled in the domain. Each cell is associated with one point in that any point in such a cell is almost closer to a point than to any other point in . The word “almost" is used to indicate the exceptions where this point may be equally close to two or more points in . Such a relationship can be described as where are the coordinates of and is the Euclidean distance between and any point at interior of .
The Voronoi diagram is directly constructed from a set of points by using the software Matlab as illustrated in Figure 1. The coordinates of these points are vital because they control the shapes as well as the area of generated polygons. In order to guarantee the areas of those arbitrarily generated polyhedrons within a prescribed range, the distance of any two spatial points in must be checked to meet One simple method to generate those sets of points is to insert randomly generated points sequentially until they eventually saturate the whole domain. When a new point is to be inserted, (6) must be checked against all the accepted points. This process continues until the domain is nearly saturated, and the probability of the acceptable point becomes low and further trial point is rejected. The Voronoi diagram can then be constructed based on these selected points.
(a) Random spatial points
(b) 2D Voronoi diagram
3.2. Artificial Components Generation
After the Voronoi diagram is ready, the generated polygons can be classified to represent different types of minerals. A series of pseudorandom numbers are generated corresponding to the generated polygons one by one. According to the statistical contents of the minerals in the polycrystalline rock, these pseudorandom numbers can be divided into several groups. Each group corresponds to one type of mineral. Thus, each polygon is specified to be one type of artificial component according to its associated pseudorandom number. Once the ratio of these generated artificial components meets the prescribed one, such a process stops. It is considered that the mechanical properties in the same mineral component are not homogeneous due to the inherent preexisting defects. A treatment with spatial distribution of the mechanical properties within the individual mineral component is still necessary. For example, the strength variations at the particles over the mineral component can be incorporated in the analysis by assuming their properties to obey a statistical Weibull distribution. Each kind of components can have an independent degree of heterogeneity. Thus, the material heterogeneity is reflected at both mineral level and particle level.
3.3. Representation of the Physical Domain Microstructure by the SPH Particles
In order to use the SPH particles to represent the physical domain’s microstructure, one must determine to which polygons each particle belongs. If the particle’s center falls inside the polygon, it is assigned to the component’s properties represented by this polygon. There are many methods to judge whether a point is inside or outside a 2D polygon, such as the “crossing test” method by Shimrat [15], the “triangle test” method by Didier [16], and the well-known “angle summation test” method. Implementations of these methods can be referred to in Haines [17]. The current work employs the “crossing test” method for its efficiency. If the number of the crossings of the ray and the polygon’s edge is even, the test point is outside the polygon. Otherwise, it is inside. Figure 2 gives an illustration of such a process by using the regularly arranged particles.
3.4. 2D Granite Microstructure Generation
The granite is generally light grey and medium grained which can be discerned easily by the naked eye. Statistical data from Geological Unit of Singapore [18] show that the Bukit Timah granite in Singapore is mainly composed of the feldspar, quartz, and biotite grains and their volumes are in a ratio of around 6 : 3 : 1. The grain size generally ranges from 3.0 mm to 5.0 mm. Table 1 presents the mineral contents and their properties in the Bukit Timah granite. Figure 3 shows one cross-section image of such a granite sample. One rectangle artificial rock specimen has been constructed based on the above-described approaches to resemble the Bukit Timah granite. The width of this specimen is 5.0 cm and the height is 10.0 cm. It is discretized by the regularly distributed 31250 SPH particles with the smoothing length of 0.4 mm.
The spatial points prepared for the Voronoi diagram are firstly generated to fill the domain according to (6). In the equation, the minimum and maximum lengths are specified to be 0.3 mm and 0.5 mm, respectively, to meet the requirement of the typical granite grain size. After these polygons are readily constructed, they are classified into three different mineral types and the heterogeneities are assumed to follow the Weibull distributions. Figure 4 shows one typical granite sample, named N1. The grains in dark black are biotite. Quartzite grains are in grey white and the feldspar ones are in grey. Mineral contents in this artificial granite sample are presented in Table 2. As can be seen, they are very close to the expected ratio in terms of the volumetric contributions.
3.5. A Unified Twin Shear Strength Criterion
An elastoplastic damage model that can effectively represent the mechanical behavior of rock-like material failure is adopted. It is comprised of two surfaces to represent the strengths of intact and fractured materials. A damage model describes the evolution of the material from the intact state to the fractured state. The strength criterion is based on an extension form of the unified twin shear strength criterion (UTSS) by Yu et al. [19], which includes two hydrostatic pressure dependent meridians representing the generalized tensile and compressive strength states, respectively. According to the UTSS criterion, when , it gives the upper bound of the convex failure surface; when , it yields the lower bound. And any failure criterion satisfying the convex requirement can be derived by a linear combination of the lower and upper bounds using the parameter . Yu et al. [19] suggested that, for rock materials, can take a value between 0.5 and 1.0. Detailed mathematical representation and the materials parameters of the UTSS criterion have been given in Ma et al. [20, 21].
4. Numerical Model
4.1. Specimen Geometry and Loading Conditions
Based on the previous generated artificial specimen N1 by the proposed microstructure modeling method, two granite samples are cut by intercepting the central part of specimen N1 as shown in Figure 5. Of them, one is used to conduct the Brazilian splitting test simulation with a diameter of 0.05 m; the other is for the uniaxial compression test with a width of 0.05 m and a height of 0.03 m. The mineral contents for these two specimens are listed in Table 2.
In both simulations, the samples are sandwiched between two rigid walls with the prescribed velocities defined as where . The velocity ramp can help to achieve the stress equilibrium and prevent the premature failures due to the stepwise increasing force on the loading boundaries. After , the prescribed velocities on the boundary are kept at a constant value . The strain rate is constant after the time which is calculated by , where is the nominal length of specimen along the loading direction. The loading cases and their corresponding strain rates for these two kinds of tests are presented in Table 3.
4.2. Numerical Approach and Material Parameters
The developed SPH code is applied to conduct the Brazilian splitting test and uniaxial compression test simulations under the plane stress condition. Model parameters follow those of their “parent” specimen N1, in which feldspar, quartz, and biotite grains are assumed to have the independent Weibull distributions with the shape values of 50, 5, and 10, respectively. The SPH particles in these two specimens have the same smoothing length of 0.4 mm. The total particle numbers in the Brazilian test and compression test samples are 12281 and 9375, respectively.
5. Simulation of the Brazilian Splitting Test
5.1. Loading Force-Displacement Curves and the Fracture Processes
The compressive loading force with the loading displacement curves for the cases with the strain rate from to is shown in Figure 6. In the figure, the dashed, dotted, and solid lines are the profiles of the upper boundary force, the bottom boundary force, and their average one, respectively. The corresponding simulated fracture processes and the final failure patterns of the specimen under these strain rates are depicted in Figure 7. The failure process is presented by a series of damage distribution frames starting from the moment when the first visible crack occurs. These frames are indicated by the white circle marks on the loading force versus curve in each case as shown in Figure 6.
| (a) Strain rate |
| (b) Strain rate |
| (c) Strain rate |
| (d) Strain rate |
| (e) Strain rate |
| (a) Strain rate |
| (b) Strain rate |
| (c) Strain rate |
| (d) Strain rate |
| (e) Strain rate |
For the two cases with the strain rates and , their loading forces versus displacement curves are similar. The loading force increases almost linearly until the peak load and then drops suddenly once the sample loses its bearing capacity. Their fracture processes are almost identical. It can be seen that the crack is initiated near their peak loads on the quartz and biotite grain boundaries around the disc center, as indicated by those open circles. Due to the effect of specimen’s heterogeneous microstructure, the initiated crack may not occur at the disc center but take place where the stress condition is most severe as discussed by Andreev [22]. With further increasing load, the crack propagates rapidly towards the two ends of the disc along the loading diameter until it cuts through the whole disc. As can be seen from the failure pattern, the predicted crack splits the disc into two halves.
For other higher strain rate cases, the shapes of the loading force curves differ remarkably from those of the two lower strain rate cases. In these higher strain rate cases, the curve exhibits somewhat nonlinear characteristic after the first visible crack occurs, as marked by the first white circle. For the case with a strain rate of , the first visible crack occurs around 87% of its peak force and then propagates towards the two ends. At its peak force, another crack appears along the boundary between quartz and biotite grains. The propagations and interactions of these two major cracks result in the oscillation of the response force in the subsequent stages. The failure processes in the cases of the strain rates and are very much similar. In both cases, the tensile crack is initiated much earlier before their peak forces (around 58% in the strain rate case and 61% for the case). Their failure processes can be featured by the developments and interactions of many short cracks near the loading diameter plane. The loading force curves behave in a much more fluctuant way due to the intensive crack activities. In addition, the bending cracks can be clearly observed along the disc edge.
From Figure 6, it can be seen that, in all cases, the reaction forces on both loading sides are almost identical before their peaks. At the postpeak stage, the specimen’s failure becomes unstable and the force oscillations can be observed especially in the higher strain rate cases.
5.2. Strain Rate Effects
The predicted fracture processes shown in the previous section manifest that their failure processes are influenced by the strain rates notably. In the lower strain rate cases with the strain rates and , specimen’s failure is roughly due to the development of one major crack initiated at the weaker zone around the disc center. In those higher strain rate cases, more short cracks are developed. Consequently, the coalescence of the short cracks forms the crack band, which finally causes the specimen to rupture.
Figure 8 shows that the specimen’s macromechanical behavior differs considerably under different strain rates. In the lower strain rate loading cases with a strain rate from to , a major crack is initiated and propagates and the disc completely loses its bearing capacity shortly after its peak. However, in those higher strain rate cases, the applied forces can still increase until the isolated short cracks coalesce to form the crack band. The interactions of those short cracks during such processes cause fluctuations of the loading force curve.
As known to all, the SPH particles can separate naturally at large deformations. This advantage makes the method easier to process the fractures and fragmentations compared with other conventional grid based methods. Figure 9 shows the predicted failure patterns in the strain rate of , , and , respectively. Those blank areas in the plotted graph are cracks formed by the particle separations. As can be seen, the failure patterns resemble those obtained in the quasistatic and SHPB dynamic experimental tests. When the strain rate approaches the quasistatic condition, the specimen is split into two pieces by the primary crack. However, at a higher strain rate, more fragments and conspicuous crashed zones are formed near the two loading boundaries.
(a) Brazilian splitting test
(b) Granite compression experimental results under the static test and SHPB tests
6. Simulation of the Uniaxial Compression Test
6.1. Axial Stress-Strain Curves and Fracture Processes
Figure 10 plots the specimen’s axial compressive stress-strain curves with the strain rate from to . In each case, the upper boundary stress, the bottom boundary stress, and their average are plotted in the dashed, dotted, and solid lines, respectively. The corresponding fracture processes and the final failure patterns of the specimen under these strain rates are depicted in Figure 11. Similarly, the first failure process frame in each case corresponds to the first visible crack in the specimen. These frames are indicated by the white circle marks on the curve shown in Figure 10.
| (a) Strain rate |
| (b) Strain rate |
| (c) Strain rate |
| (d) Strain rate |
| (e) Strain rate |
| (a) Strain rate |
| (b) Strain rate |
| (c) Strain rate |
| (d) Strain rate |
| (e) Strain rate |
For cases with the strain rate from to , their stress-strain curves are similar resembling a typical brittle failure property. The axial stress increases almost linearly until its peak value and then drops suddenly. Their failure processes are almost identical. The crack is initiated at the peak stress at the boundary between the quartz grain and biotite grain, as indicated by the open circle. The crack propagates rapidly at the postpeak stage and coalesces with another developed crack. Finally, these cracks form one vertical macrotransgranular crack across the specimen. More cracks occur and propagate rapidly in the direction parallel to the loading axis. Besides, kinking in suitably oriented biotite grains can be evidently observed. Some of them form the small shear faults as shown in the specimen’s failure pattern plot. For other higher strain rate cases, their axial stress-strain curves exhibit three recognized regions: the linear elastic, the nonlinear, and the postpeak softening region. In the case of strain rate , the initiated crack occurs around the 98% of its peak stress. The macrocrack propagates across the specimen at 98% of its peak stress at the postpeak stage. For the strain rates and cases, their fracture processes resemble one another. In both cases, the cracks appear much earlier compared with other lower strain rate cases (around 85% of its peak stress in the strain rate and 75% in the cases). The major cracks are completely formed at their peak stresses. In the subsequent postpeak stage, more short cracks are developed predominantly along the loading direction. These cracks, together with the kinking in biotite and other grains, form several faults and divide the specimen into many fragments.
6.2. Strain Rate Effects
The compressive axial stress-strain curves under different strain rate cases are plotted together in Figure 12, where the stress is taken as the average between the upper’s and bottom’s. Result shows that the specimen’s macromechanical behaviors are remarkably affected by the strain rate. It is clearly observed that the maximum stress and the maximum sustained strain of the specimen increase when the strain rate is getting higher. When the strain rate is lower, the sample experiences little plasticity before a brittle failure occurs. The stress reaches the peak linearly and drops down sharply after its peak, whereas, for the higher strain rate cases, the curves experience linear elastic, nonlinear, and the postpeak softening stages.
Figure 13 shows the simulated final failure patterns of cases with strain rates , , and , respectively. As can be seen, the sample fractures into more small fragments when subjected to a higher strain rate loading.
7. Strain Rate Effect on the Predicted Dynamic Tensile and Compressive Strengths
The predicted dynamic tensile strength () and compressive strength () under the different strain rates are presented in Table 4. is calculated using (8) by taking the maximum average compressive force on two loading boundaries as the peak force. directly takes the maximum average response stress on the two loading boundaries, As can be seen, both and exhibit the strain rate sensitivity in a sense that they increase as the applied strain rate increases. Such a relationship is coincident with the results reported by many researches. However, to quantitatively validate these predicted results with the experimental results is difficult. It is simply because the artificial specimens used in the current analyses are constructed by the statistical heterogeneous microstructure model. Although their microstructures have been entitled with the major features of granite samples, these artificial specimens are not the replicas of those samples used in experimental tests. In addition, the mechanical parameters of the mineral components are based on an approximation method with some assumptions due to the difficulties in obtaining them from measurements directly.
Figure 14 compares the numerically predicted dynamic tensile and compressive strengths with those obtained from the experimental tests on the granite samples, respectively. The comparisons take the forms of their dynamic increment factor () for the strengths. The value is calculated by normalizing the dynamic strength with respect to their quasistatic values (at a strain rate around ). Due to the expensive computational effort to perform the quasistatic analysis using the SPH code, the values of the numerically predicted strengths are normalized with that at the strain rate case instead. Similarly, the values of the specimen’s average pressure () in the numerical simulations are also calculated under different strain rates. As seen from the two plotted diagrams, the experimental and numerical results demonstrate similar trends in the sense that both the tensile and compressive strengths increase with the applied strain rate. It is also observed that the increments of numerically predicted strengths are generally not that sharp as those obtained by experiments. Such deviations may come from several sources. One is due to the limitations in the current model as explained previously. Another thing is that the values of the predicted strengths are normalized by the higher strain rate values instead of those at the quasistatic conditions. This will probably underestimate the predicted values.
(a) Tensile strength
(b) Compressive strength
It is very interesting to observe that the trends of of tensile and compressive strengths curves remarkably match those values of their average pressure, as shown in Figure 14. Such a good agreement clearly attests that the pressure increment caused by the strain rate plays an important role for the rate-dependent dynamic strength. Bedsides the lateral inertial effect, the crack activities can also be one of the factors to cause the increment of the specimen’s internal pressure. As depicted in the fracture processes of the Brazilian splitting and uniaxial compression simulations (see Figures 7 and 9), the specimen’s failure is closely related to the activities of the cracks, including their initiations, propagations, and coalescence. At the lower loading rate cases, only few microcracks are activated before the maximum load (or stress) is reached. Finally, one single macrocrack causes the specimen a catastrophic failure. On the contrary, at the higher loading rate cases, it is clear that substantially more microcracks are developed before the maximum load. Because the extension of the microcracks takes time, other microcracks can be activated before a neighboring activated microcrack extends and unloads them. Consequently, the specimen breaks into more pieces. The formation of these cracks will not only delay the specimen’s failure but also cause the stress release and consume considerable energy.
8. Conclusions
In this paper, a microstructure modeling method is proposed to consider the spatial variation of the components distributions as well as their strength distributions of multiphase material. This method has been successfully applied to construct the artificial granite specimen and the granite rock dynamic failures are further examined by using the artificial specimen. Numerical simulations on Brazilian splitting tests and uniaxial compression tests have been conducted using the developed SPH code. The numerical results are compared with those from the experiments in terms of the final failure patterns and the dynamic strength increments. Results demonstrate that they have good qualitative agreements.
The predicted fracture sequences in both Brazilian splitting tests and uniaxial compression tests manifest that their failure processes are notably influenced by the applied strain rates. In the lower strain rate cases, the failure is roughly due to the development of few major cracks initiated at a weaker zone. On the contrary, under the higher strain rate cases, the specimen’s failure resulted from the coalescences and interactions of many short cracks.
Specimen’s macromechanical behaviors exhibit pronounced strain-rate sensitivities. As the strain rate becomes higher, both the specimen’s loading capacity and its sustained maximum displacement become larger. Finally, the strain rate effect on the strength is analyzed and compared with the experimental results. It seems that such an effect is attributed to the internal pressure changes due to the inertial effects as well as the specimen’s crack activities under different strain rates.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.