Abstract
As a Lagrangian mesh-free numerical method, the Smoothed Particle Hydrodynamics (SPH) method has been traditionally applied for modeling astrophysics, fluid flows and thermal problems, and there has been a growing interest in applying SPH to solid deformation problems. However, the potential of this method for quasistatic analysis of rock-like brittle materials has not been clearly explored. The major aim of this paper is to investigate the effects of key factors in SPH on the load-deformation response of rock-like solids, including variations in the particle approximation theory, the magnitude of the smoothing length and its variable method. Simple uniaxial compression (UC) loading conditions were chosen, and a series of numerical studies were carried out sequentially on an idealized elastic case and an actual test of marble material. Typical results of the axial stress-strain response from infinitesimal to finite deformation as well as the progressive failure process for the marble tests are given and the influences of various factors are discussed. It is found that only provided proper choices of particle momentum equation and the smoothing length parameter, the SPH method is capable for favorably reproducing the deformation and progressive failure evolution in rock-like materials under quasistatic compression loads.
1. Introduction
SPH is a mesh-free, adaptive, Lagrangian particle method that can be used to obtain solutions to systems of partial differential equations. In contrast to the concept of discretization methods which discretize a continuum into a finite set of nodal points, SPH consolidates a set of discrete particles into a quasicontinuum, and each particle represents specific material volumes. Since there is no mesh to distort, the method can handle large deformations in a pure Lagrangian frame [1], hence presenting extraordinary power for easy treatment of void regions, material interfaces, and multifracturing in materials. The technique was initially formulated to solve astrophysical problems [2, 3] and has been widely adopted in fluid dynamics related areas [4–6]. Material strength can be considered in a fairly straight forward manner, and Libersky and Petschek [7] were the first to incorporate an elastic-perfectly plastic material strength model into the SPH framework and they applied it to the simulation of the impact of an iron rod with a rigid surface. In recent years, there has been a growing interest in applying SPH for modeling solid deformation problems in a broad range of applications. However, most of these studies were focused on the high velocity impact, blasting, and granular flow types of problems, such as the study of the deformation of a metal cylinder resulting from the normal impact against a rigid surface [8, 9], modeling of impact induced fractures in brittle solids [10, 11], and the simulation of broken-ice fields floating on the water surface and moving under the effect of wind forces [12]. Rather limited publications can be found to date with regard to the application of the SPH method in quasistatic analysis of deformation and possible failure response of solid materials.
Given the extraordinary potential of SPH in dealing with large deformation, the method has also obtained a number of applications in the field of geotechnical engineering. Previous applications of the SPH method in geotechnics mostly dealt with fluid-like flow problems, in particular, those involved in analyses of dam-break flood hazard [13, 14], water-structure interaction [15, 16], and pore-scale flow phenomena in soil-like porous media [17, 18]. More recently, Bui et al. [19] also developed SPH to solve the large deformation and postfailure flow of noncohesive and cohesive soils. Ma et al. [20, 21] applied SPH to investigate the static or dynamic failure of brittle heterogeneous materials by tracing the propagation of the microscopic cracks as well as the macromechanical behaviors, in which the material heterogeneity was modeled by a statistical approach.
Furthermore, with the fast development of SPH theory, quite a number of extensions using different particle approximation theories, varied choices of the smoothing length and correction techniques have been proposed. To the best of our knowledge, no discussion of various options in SPH and their effects on the modeling of progressive deformation and in turn failure process of rock-like materials under quasistatic loads has been reported in literature. The primary aims of this paper are to study the suitability of the SPH method for quasistatic modeling of rock-like materials and to study the effects of influencing factors on the load-deformation response and the associated failure progress. The paper is organized as follows: firstly, a brief introduction to the fundamentals of SPH is provided, forming the basis of following parametric numerical studies and discussions; secondly, the main part of this paper presents a series of numerical investigation into the effects of varied options in SPH on the uniaxial compression response, in which two typical cases are chosen, an idealized elastic case of a cubic specimen and an inelastic case considering possible damage and failure of marble materials; in the final part, a summary of the parametric studies and some conclusive discussions are provided.
2. Fundamentals of SPH Method
2.1. Basic Formulations
The theoretical fundamentals of SPH can be described in two main steps. The first is the integral representation or kernel approximation of the field functions. The second step is the approximation of particle variables, that is, the modeling domain is discretized into a finite number of particles which carry field variables such as mass, density, stress, and so forth and move along the particle velocity. As in the framework of continuum mechanics, the governing equations for constructing SPH formulations are mainly composed of mass, momentum, and energy conservation laws.
The material properties of each particle, such as velocity, density, stress, and so forth, are calculated through the use of an interpolation process over its neighboring particles, which is based on the integral representation of a field function as follows: where denotes the kernel or smoothing function, and is the smoothing length, which defines the influence domain of . Generally, the chosen kernel function should satisfy three conditions: the normalization condition, the delta function property, and the compact support condition as follows: in which is a constant that specifies the nonzero region of the smoothing function for a point at a position vector . This implies that the integration over the entire domain is localized to an integration over the support domain of the smoothing function.
The integral representation (1) can be then discretized as a summation over the particles within the support domain as follows: where denote particles within the support domain of the particle at , , and are the mass and density of particle , respectively. In a similar manner, the approximation for the spatial derivatives can be obtained.
Using the above particle approximations for a function and its gradient at a particle, the conservation laws in the form of partial differential equations can be then converted into equations of motion of discrete particles and then solved by an updated Lagrangian numerical scheme. However, different transformations or operations may result in different discretized forms of SPH equations [19, 22]. In the coming section, some typical and commonly used SPH formulations shall be presented, forming the basis of subsequent discussions on the numerical results of the chosen cases. More details on the SPH method and technical treatments can be found in the textbook by G. R. Liu and M. B. Liu [23] and a recent review of SPH development by Monaghan [24].
2.2. Choices of Particle Approximation Theory
A large variety of SPH formulations have been proposed by many researchers, which originate from the application of different forms of kernel functions and variable transformation options. Obviously the choice of the kernel function in (1) will directly influence the accuracy, efficiency, and stability of the numerical analysis [22, 23]. Several types of kernel functions have been proposed in literature, such as bell-shaped function, Gaussian function, quadratic function, and spline functions of cubic or higher order. By far the most widely used kernel function in SPH is the one devised by Monaghan and Lattanzio [25] based on the cubic B-spline function, which is given as for three dimensions where is the normalized distance, and the coefficients assure proper normalization. This function has the advantage of having compact support and a continuous second derivative and was chosen for the following numerical analyses.
As the main focus of this paper is on quasistatic deformation and progressive failure analysis of rock-like solids using the SPH method, obviously the equation of motion plays a dominant role during this type of analysis. Hence, our attention is given to the variable choices of discrete interpolations for the linear momentum equation. Four types of typical formulations are chosen in this study and outlined as follows. It should be noted that the formulations given in this section are all constructed assuming no body force, no mass, and no heat sources involved.
Formula I:
Formula II:
Formula III:
Formula IV: where denotes the velocity of each SPH particle, and the indices and denote the particle number; is the stress tensor; denotes the gradient of kernel function , which is taken with respect to the as . and , respectively, denote the gradient of kernel function with respect to the and . One can easily note that Formulas III and IV are established through a renormalization correction on the former two equations, respectively, as denoted by a second rank tensor in (7) and (8). The purpose of the correction treatment is to impose general boundary conditions as well as to diminish the effect of particle deficiency at boundaries due to inaccurate sum estimates. Detailed formulation and technical treatments of the tensor can be found in [1], which extended the “ghost particles” approach by Libersky and Petschek [7].
The above Formula I represents a general summation interpolant for the linear momentum equation. Specifically with an assumption that the smoothing function is symmetric and the smoothing length at each particle remains the same constant, the equation will be equivalent to Formula II. Also some other summation interpolants can be found in literature using different approximation rules.
Commonly the continuity equation and the energy equation are also transformed into summation interpolants as parts of fundamental SPH equations in continuum mechanics. Using the same approximation rules as in the derivation of the above momentum equation, the SPH continuity and energy equations can be given as follows: where is specific internal energy, denotes the hydrostatic pressure for the stress tensor at a particle. It should be noted that the use of continuity and energy equation may be optional in the current SPH application to quasistatic deformation and failure analysis of rock-like materials.
2.3. Choices of Smoothing Length
The smoothing length represents the influence width of the kernel and its magnitude determines the number of particles with which a given SPH particle interacts. Hence it is straightforward to raise a problem about the influence of the value on quasistatic modeling of deformation and failure process of rock-like materials. During the explicit integration of the above governing equations, routinely the initial value of the smoothing length for each particle is determined as the maximum of the minimum distance between every particle. Moreover, a scaling factor denoted by is commonly introduced for the adjustment of value, which means that the actual influence width of the kernel will be scaled to . Obviously a choice of value should be not less than 1.
Another important aspect in the choice of is about its evolution during the entire simulation. Early SPH simulations used a fixed smoothing length for all particles, which indicates that the number of particles within the influence zone significantly depends on the type of loading. For instance, for compressive loading condition the number of particles increases, whilst for tensile loading less particles would be inside the influence zone. However, the above summation interpolants imply that the accuracy of an SPH simulation depends on having a sufficient number of particles within the smoothing length. Allowing each particle to have its own variable smoothing length which changes with time and space according to local conditions may increase the spatial resolution substantially [26, 27]. Different ways for the evaluation of the variable smoothing length have been proposed by many researchers. Two types of frequently used definitions of the variable are given as follows.
method I:
method II: where means the divergence of velocity . The smoothing length increases when particles separate from each other and reduces when the concentration of particles is important, aiming to keep the same number of particles in the neighborhood.
3. Numerical Experiments of UC Test
In this section a simple unconfined compression test was chosen and parametric studies were conducted. The numerical investigation was comprised of two parts as follows: in the first part, a series of elastic analyses were conducted on an idealized cubic sample as a benchmark case, through which the influences of various SPH options on the load-deformation response from infinitesimal to large axial strain conditions, were analyzed and discussed; in the second part, inelastic analyses of typical UC tests of Georgia Cherokee marble material by Hudson et al. [28] were carried out considering possible damage and failure behavior. Three-dimensional numerical studies were chosen for these two types of analyses using the commercial software LS-DYNA ver971 [29, 30], in which the above-described particle approximation formulations and the control of magnitude as well as its variable methods have been well implemented.
3.1. Idealized Elastic Case
3.1.1. Model Description
For the idealized elastic case, a cubic specimen of 100 mm in side length was chosen as a threefold symmetric distribution of SPH particle can be obtained along three orthogonal directions. As shown in Figure 1, the numerical models consist of three major components, namely, (1) the upper rigid platen, (2) the lower rigid platen, and (3) the specimen in between. The dimensions of the two platens were defined to be a bit larger in section (120 mm × 120 mm) and 10 mm in thickness, such that a full contact between the platen and the specimen end surfaces can be maintained even at the stage of large lateral deformation induced by top compression.
| (a) Type-A ( particles) |
| (b) Type-B ( particles) |
| (c) Type-C ( particles) |
In the numerical model, the upper and lower steel platens were modeled by linear hexahedron solid finite elements (FEs). The Young’s modulus of the platen was taken as MPa, and Poisson’s ratio equal to 0.28 as commonly applied values for steel material. For the cubic specimen modeled by SPH particles, the elastic constants were chosen as: Young’s modulus = MPa, and Poisson’s ratio = 0.2, similar to those of common types of rock material. Regarding the mass properties of the platens and the specimen, the magnitude of density parameter would not influence the load-deformation response for such a quasistatic analysis. However, the choice can greatly affect the allowable time increment and in turn the computational cost. Herein the density of the platens was defined as kg/m3, a common value for steel material. For the specimen, the density was taken as kg/m3 as common types of rock material, and a mass scaling of was applied. The Courant-Friedrichs-Lewy (CFL) condition is necessary for convergence in the explicit integrations of the above system equations, which requires the time step to be negatively proportional to the sound of speed. Obviously the above mass scaling can enlarge the allowable time increment by about 1000 during the explicit integration analysis.
Regarding the load and boundary conditions, the bottom platen was specified as a fixed support. The compression was defined by a vertical displacement control on the top platen. A final magnitude of 10 mm was chosen for obtaining a large axial strain level of about 10% finally, such that the suitability of the SPH method for modeling elastic solids from infinitesimal to large strain can be examined. The top compression was initiated from zero and increased linearly to 10 mm in a time period of 30 seconds, which was sufficiently slow for the requirement of a quasistatic analysis through a posterior check. For the interaction between the SPH particles and the FE domain, a master-slave contact algorithm was applied to couple the two techniques. In this study, a penalty based “nodes to surface” contact [29] is adopted. The SPH particle elements are defined as the slave side and the finite elements are defined as the master side. Zero friction contact interaction was defined between the specimen and the two platens, simulating an ideal unconfined compression case.
3.1.2. Analysis Programme
Four series of analyses have been conducted and the details of each series are listed in Table 1. Series 1 aims to examine the effect of particle density, using three types of particle distribution models as given in Figure 1. The effect of different particle momentum equations, as given in the above section, is examined in Series 2 using the Type-B particle distribution model. Series 3 investigates the influence of the variation in magnitude by adjusting the parameter, and three values falling in a range of , as commonly adopted in literature, are considered. The analyses in Series 4 aim to study the effects of the two variable methods during the integration process, which are, respectively, shown in (10) and (11).
3.1.3. Results
A summary of the axial load-displacement results from the series of SPH analyses is given in Figure 2, and the individual effect of the above factors on the UC response for the cubic specimen is demonstrated. Moreover, a three-dimensional single element FE analysis for the UC test has been conducted and the modeling result is also shown in Figure 2 for comparison. It can be observed that generally all the predictions using SPH method with varied options in Table 1 compare favorably with the one-element FE solution at an early stage with relatively small top compression, lower than 1 mm for the studied model, whilst more noticeable discrepancy is shown at a later stage with greater compressions applied.
(a) Effect of particle density
(b) Effect of particle momentum equation
| (c) Effect of parameter |
| (d) Effect of variable method |
The effect of SPH particle density is shown in Figure 2(a). It is found that at the early stage with low value of top compression, approximately at a level of axial strain lower than 1%, the predictions by the three types of models are close to the analytical solution, all showing a linear increase of axial load with respect to the top compression. A comparison of the slope of the curves at this stage shows that amongst the three models, the Type-C model using the finest particle distribution gives the best prediction with a minor difference (~4.7%) with respect to the one-element FE solution, and the predictions by the other two models using coarser particle distribution give a slope difference of 8.6% (Type-A) and 6.1% (Type-B), respectively, relative to the FE solution. One can also note that with a greater than twofold increase of SPH particles from 27000 (Type-B) to 64000 (Type-C), only a small improvement in the accuracy of load-displacement response is shown by the SPH method.
With a larger top compression applied at the later stage, geometric nonlinearity plays a more significant role and the stiffening effect induces a slightly concave-downward curve. It can be seen that all the modeling results in Figure 2(a) reproduce this essential behavior, also one can note that comparatively more prominent discrepancy is shown between the simulation results and the FE solution with the axial deformation increased, particularly at an axial strain level greater than 2.5%. Again, it is shown that the Type-C model with the finest particle distribution gives the best prediction at a larger top compression, except at an intermediate stage with the top compression falling in the range of 1 mm to 2.5 mm due to the occurrence of a nearly horizontal plateau in the curve. The underlying mechanism of this phenomenon is addressed in the Discussions section.
Figure 2(b) illustrates the influence of varied choices of the particle momentum equation. It is shown by the comparison that a better prediction can be obtained by Formulas III and IV in which the renormalization technique was applied, and the two choices basically gave the same results. For the other two choices using Formulas I and II, the simulation results are basically the same and compare less favorably with the one-element FE solution. Also it can be observed that an occurrence of nearly horizontal plateau is shown by all the SPH modeling curves at a compression level of around 1 mm.
The effect of variation in by adjusting the parameter is shown in Figure 2(c). It is found that with the parameter increased from a standard value 1.0 to 1.2, the predicted load-displacement curves were significantly influenced at the later stage with a larger top compression. In particular, for the case with , a top compression greater than 4.5 mm could trigger a significant deviation from the exact solution, which seems unfavorable and the mechanism behind this is addressed in the Discussions section.
A comparison of the results using the two types of variable methods during the integration process is shown in Figure 2(d). The two curves are approximately the same from the initial infinitesimal compression to the final large deformation, which illustrates that the two types of variable methods in (10) and (11) may induce only negligible influence on the load-displacement response for the UC problem in this study.
3.2. Inelastic Case of Marble UC Tests
3.2.1. Model Description
To further investigate the capability of the SPH method as well as the effects of the above factors on deformation and failure analyses of rock-like materials, parametric studies were conducted in this section for typical UC tests of Georgia Cherokee marble material by Hudson et al. [28]. In their study, quite a number of UC tests considering a wide variety of size and shape have been conducted for the specified marble material. The cylindrical test samples with a diameter of 101.6 mm (4 inches) and a height/diameter ratio equal to one were chosen in this study. Referring to the above elastic analyses, a particle discretization of 30 lines of particles along the three orthogonal directions, the same as that in Type-B model in Figure 1, was considered fine enough for the following analyses. The established numerical model was shown in Figure 3. A total of 21480 SPH particles were defined to approximate the cylinder specimen. Also one can note a regular distribution along the horizontal section was adopted rather than a circumferentially equally spaced distribution, which was applied to minimize the inter-particle distance discrepancy and better approximation can be expected [29]. The upper and lower steel platens were modeled by finite elements and their diameters were defined the same as that of the specimen in accordance with the actual test conditions [28]. A constant friction coefficient was applied for the definition of the contact interaction between the platens and the marble specimen, as no special measure for diminishing the frictional restraint effect was applied during the tests.
(a) Side view
(b) Horizontal section view
A coupled damage-plastic material model by Malvar et al. [31] was adopted to simulate the marble behavior. The model supports nonlinear elasticity and uses a three-invariant formulation for the failure surface. It is characterized by the employment of three failure surfaces, including the maximum, residual, and initial yield surfaces, for a proper representation of material behavior along multiple stress paths, including uniaxial, biaxial, and triaxial tension and compression. Although the model was originally developed for modeling the deformation and damage behavior of concrete materials, it was chosen in this study based on the presumption that the constitutive response characteristics of most types of rock and concrete materials are relatively similar. More details about the constitutive model can be found in [31].
Regarding the inputs of relevant parameters for the Georgia marble, as limited test data of the UC axial stress versus strain curves were provided in [28], only the deformation modulus and the UC strength parameters could be calibrated from the test results. Even though a variety of numerical approaches have been proposed by many researchers for the parameter calibration based on incomplete test data, such as those based on artificial intelligence [32], the calibration process in this study was simply based on trial analyses. The initial values of most parameters for the above model were determined using the internal parameter generation capability by Malvar et al. [33], during which only the unconfined compressive strength is needed and typical concrete test data are taken as a solid basis. Some minor adjustments were further made on the autogenerated parameter values as follows: a uniaxial tensile strength of 4.6 MPa for the same type of marble given by Schwartz [34] was adopted; and a slight adjustment of the compressive scaling exponent ( parameter) from the autogenerated value of 1.6 to 1.0 was applied for a better match of the observed UC softening response. A list of the input parameter values for the above material model was given in Table 2. Note that for the definition of nonlinear elasticity, the internally generated relation curve of loading and unload/reload bulk modulus with respect to volumetric strain was directly used as given in Figure 4. The corresponding shear modulus can be calculated by the bulk modulus curve and the specified Poisson’s ratio according to the prescribed procedure in [31].
3.2.2. Analysis Programme
With the analysis results of the above elastic case forming a basis, two series of analyses have been conducted for this inelastic case: firstly, the effect of the two types of particle momentum equations with renormalization technique (Formulas III and IV) on the inelastic load-deformation response was examined, as either choice gave more favorable predictions and present only slight difference on the load-deformation response in the above elastic case; secondly, the influence of the two variable methods during the integration process, as shown in (10) and (11), was investigated. Regarding the magnitude, a constant value of 1.0 was deemed most appropriate from the above elastic analysis results and was chosen in all the simulations of this section.
3.2.3. Results
A comparison of the nominal axial stress versus strain response from the numerical results is given in Figure 5, and the effects of the above two factors are separately illustrated. The corresponding test results from different sized specimens but with the same ratio equal to unity are shown in the figure. The experiments exhibit an initial nonlinearity due to the closure of microcracking within the specimens that has not been included in the models, and removed in the comparison shown in Figure 5. It is shown that generally the SPH based modeling can adequately reproduce the complete axial stress versus strain response. The alternative choice of the prescribed particle momentum equations and the variable methods would not essentially change the modeling results, and only small difference was shown at the postpeak stage.
(a) Effect of different particle momentum equations
(b) Effect of different h variable methods
All the four simulation cases gave approximately the same peak stresses, 84.4 MPa, which closely match the test results (90.7 MPa) of cubic specimens with a diameter of 101.6 mm. However, one may note the difference between the predicted peak stress and the input value of unconfined compressive strength for the marble material, 70.7 MPa given in Table 2. It should be pointed out that the input UC strength was chosen from the test data of samples with a largest ratio of 3 : 1 and a largest diameter of 101.6 mm, as routinely recommended in the suggested methods for the measurement of UC strength of rock materials [35]. It is also found that the test data by Hudson et al. [28] present an obvious trend of increase of the UC strength with the decrease of the ratio from 3 to 1. Hence the observed difference between the predicted maximum stress and the input UC strength parameter can be attributed to the specimen shape effect on the compressive strength, which was caused by severely nonuniform stress distribution in the specimens with small ratios [28].
The softening response is highly affected by the mechanism of progressive structural breakdown of the test specimen [28], and the development of damage and failure within the specimen is worthy of discussion. As the modeling results of the postpeak softening response in Figure 5 for all the cases compare favorably with the test records, the case using Formula III and variable method I was taken as an instance. The progressive damage and failure evolution at typical stages along a vertical diametrical section of the specimen is given in Figure 6, as indicated by the distribution of a scaled damage measure . For the chosen coupled damage-plasticity model [31], a value increasing from 1 to 2 means the current failure surface is contracting from the maximum to the residual failure surface (softening). As the aim of this investigation is to identify the evolution of specimen failure, a range of from 1.8 to 2.0 is selected and its distribution is shown in Figure 6. For surficial parts with approximately zero confining pressure stress, a value close to 2 would indicate the region is nearly completely damaged since the residual strength for rock material in tension is zero.
| (a) mm |
| (b) mm |
| (c) mm |
| (d) mm |
| (e) mm |
| (f) mm |
| (g) mm |
| (h) mm |
It is shown that for the marble specimen under unconfined compression, the damage initiated circumferentially from the surficial parts (Figures 6(a) and 6(b)), including the portions near the edge corner and the middle region, where less lateral confinement would come from the upper and lower platens. With the top compression progressively applied, the surficial damage zones coalesced into one circumferential part and a type of slabbing failure was formed as expected (Figures 6(c) and 6(d)). At later softening stages, the compression induced damage propagated towards the internal central portion (Figures 6(e) and 6(f)) and connected at the center of the specimen at a top compression of 0.31 mm. Also it is found that the undamaged zones close to the upper and lower platens at this stage are of similar shapes as those of the horizontally confined zone for a cubic UC specimen due to frictional constraint given by van Vliet and van Mier [36]. Further top compression induced the vertical extension of severe damage zone (Figure 6(g)), and the load capacity of the specimen dropped to nearly zero at a top compression of about 0.71 mm for the interconnected fully damaged zone caused the structural collapse of the test specimen (Figure 6(h)). Moreover, a typical distribution of damage and cracking at an advanced state of failure from the marble test results is given in Figure 7. A comparison of the simulation and test results demonstrates that the apparent slabbing failure pattern and the behind mechanism of progressive structural breakdown process for the marble UC tests can be adequately and favorably reproduced by the SPH method.
4. Discussions
4.1. Mechanism of Plateaus in Load-Deformation Curves
From the above UC modeling results of idealized elastic case, it has already been noticed that each load-displacement curve in Figure 2 presents a nearly horizontal small plateau at a varied deformation level, and the plateau would be mobilized earlier at a lower magnitude of axial deformation when a relatively denser particle distribution was applied. Also the results using a variable parameter in Figure 2(c) shows that when the factor is enlarged to a certain extent, herein as an instance, the simulated load-deformation curve deviates from the analytical solution to a remarkable extent at a top compression greater than 4.5 mm. An investigation into the particle approximation process during the UC analysis has been conducted for the cause of the mentioned discrepancies. Taking the cases in Figure 2(c) as an instance, Figure 8 presents the variation of the smoothing length for typical particles with respect to the top compression during the UC analysis. As noted by the inset figure, five particles at typical positions on the top surface of the specimen, which are in direct contact with the loading platen, were chosen considering its double symmetry.
| (a) |
| (b) |
| (c) |
Generally the results in Figure 8 for the three cases reveal that there exists a good correlation between the abrupt change of value and the occurrences of horizontal plateaus or noticeable deviation from the accurate solution in the load-compression response. With a relatively small input of value, 1.0 and 1.1 in this study, a general decreasing trend is shown by the variation of with the top compression (Figures 8(a) and 8(b)), which is reasonably attributed to the UC induced contractive response as a whole. Differently, the results for the case with (Figure 8(c)) indicate that an opposite increasing trend of value is mobilized at two particles near the centre of top surface, exactly at the same stage of compression as that of the deviation point in the load-compression response. The observation can be explained by the relatively more disordered distribution of particles at this stage as compared to the other two cases, which also implied the adverse effect of enlarging the influence domain by a larger value on the accuracy of particle approximation of the SPH method (3) for the compression analysis of rock-like solids. It is also found that the modeling results of the development using the alternative choice of the variable method are basically the same as described above (Figure 8).
Overall, it is found that the occurrences of the small plateaus or even significant deviation from the accurate response were mainly attributed to the abrupt change of at the corresponding compression stage, which in turn greatly affected the number of particles within the support domain and the approximation accuracy.
4.2. Factors Influencing Load Capacity and Failure Patterns
It can be observed that the predictions of axial load versus displacement response in the idealized elastic UC test (Figure 2), even with noticeable improvement when considering the renormalization technique for the kernel function, still underestimate the load capacity, particularly when a larger axial compression is applied. The observation can be attributed to the particle approximations in SPH (3), which, as already noted by many researchers, do not have the property of strict interpolants such that in general they are not precisely equal to the particle value of the dependent variable. Hence the essential boundary conditions cannot be treated in a rigorous way. Moreover, when the particles get relatively more disordered due to large top compression, comparatively more prominent error or discrepancy from the actual value would be expected for the simulation results.
When the possible damage in UC specimens is considered, the above simulation results for the marble tests present basically the same damage evolution processes, and the correspondence of the numerical response to the physical response by Hudson et al. [28] is good. The observations indicate that the two types of momentum equations are capable of predicting the progressive structural breakdown of the UC specimen when the renormalization technique is applied. However, through additional analyses considering the varieties in the value, it is found that the damage and failure evolution can be strongly dependent on the input value. Figure 9 shows the damage distribution at a typical postpeak stage predicted by an enlarged value of 1.1 and 1.2, respectively. By comparing with the result using a value of 1.0 shown in Figure 6(f), it is clearly demonstrated that significant difference in the damage pattern can be induced by increasing the value, and the case using the original influence width of the kernel seems more capable of reproducing the test response. The observation can be attributed to the nonlocal characters in the kernel approximation (3).
| (a) |
| (b) |
Even though the SPH method suffers from reduced accuracy when large axial strain is considered in the modeling of idealized elastic UC tests, the above simulation results for the marble UC tests demonstrate the capability of this method for common compression test conditions of rock-like materials. Herein with the marble UC test as an instance, it seems that a maximum strain level lower than 0.01 may be acceptable for a good estimate of the deformation and progressive failure process using the two types of momentum equations with renormalization technique incorporated, and the two types of variable methods are both appropriate for simulating such a type of quasistatic UC problem. It is also worthy to note that a proper definition of the postpeak softening properties of rock materials is of great significance in this type of analysis and should be considered appropriately according to the particle distribution.
5. Conclusions
A series of parametric numerical studies have been carried out in this study to shed some light on the applicability of the SPH method for quasistatic deformation and progressive failure analysis of rock-like materials. The numerical experiments, including an idealized elastic case and an inelastic modeling of marble UC tests, reveal that the SPH method can be a reliable tool for adequately accurate and stable simulation of quasistatic load-deformation response, provided that suitable particle approximation rule and choice of the smoothing length are defined in the SPH model. From the analysis results using the cubic B-spline kernel function, it is recommended to adopt the particle momentum equation with renormalization technique incorporated, and a radius of influence equal to , that is, , for simulating the deformation and failure process of rock-like solids under quasistatic loads. Either choice of the two variable methods for the smoothing length aforementioned in this paper can lead to acceptable results for the UC analysis. Also it is found that the occurrences of small horizontal plateaus in the load-deformation response or even prominent deviation from the accurate solution is caused by the abrupt change of value during the explicit integration process. The optional increase of the influence width of the kernel function by raising the value can lead to an inappropriate prediction of the failure evolution in rock-like solids. It is also worthy to note that the effects of these factors on the quasistatic analysis of other types of solid problems, such as tension or shear dominant cases, are still open and need to be investigated.
Conflict of Interests
This statement is to declare that, X. W. Tang, Y. D. Zhou, and Y. L. Liu, the authors of this paper do not possess any financial relationships that might bias their work. They hereby declare that no conflict of interests exists in this work.
Acknowledgments
The authors acknowledge the support by the National Natural Science Foundation of China under Grant no. 11172321 and 51109083, as well as partial support from Fundamental Research Funds for the Central Universities under Grant no. 2012ZM0091.