Application_of_Dimensonal_Analysis_in_Calibration_of_a_Discrete_Element_Model_for_Rock

Rock Mech. Rock Engng. (2007) 40 (2), 193–211 DOI 10.1007/s00603-006-0095-6 Printed in The Netherlands Application of Dimensional Analysis in Calibration of a Discrete Element Model for Rock Deformation and Fracture By A. Fakhimi1;2 and T. Villegas3 1 Departments of Mineral and Mechanical Engineering, New Mexico Tech, Socorro, NM, U.S.A. 2 Department of Civil Engineering, University of Tarbiat Modarres, Tehran, Iran 3 Department of Civil and Mining Engineering, University of Sonora, Hermosillo, Sonora, Mexico Received January 15, 2005; accepted March 18, 2006 Published online June 9, 2006 # Springer-Verlag 2006 Summary A discrete element approach was used in the simulation of rock fracture. The numerical synthetic material was made of rigid circular particles or cylinders that have interaction through normal and shear springs. The cylinders were bonded to each other at the contact points to withstand the applied loads. To characterize the microscopic properties of this synthetic material, a dimensional analysis approach was presented. It was shown that the dimensionless parameters and graphs obtained were useful tools for fast and efficient calibration of a synthetic material. This calibration method was employed for finding a numerical model for Pennsylvania Blue Sandstone. The numerical model could mimic many deformational and failure characteristics of the sandstone in both conventional and some non-conventional stress paths. Keywords: Rock fracture, discrete element method, dimensional analysis, calibration. List of Symbols D specimen width and specimen diameter E Young’s modulus E0 apparent Young’s modulus kn normal stiffness of a contact ks shear stiffness of a contact nb normal bond of a contact p maximum applied load in a Brazilian test and average of maximum and minimum principal stresses q half of difference between maximum and minimum principal stresses 1. Introduction The issue of rock deformation and fracture has been studied intensively by researchers in the past. Two general paths are usually followed to investigate rock fracture, a continuum model and a discrete model. In the continuum model, different techniques are implemented; the most common ones are softening plasticity and damage me- chanics models. In the softening approach, rock friction angle or cohesion can be considered to evolve by an invariant of plastic strain or plastic work (Vermeer and De Borst, 1984), while in the damage model, rock stiffness is reduced with deformation (De Borst, 2002). In the continuum framework, the geomaterial failure and localiza- tion of deformation can be studied as a bifurcation phenomenon as well (Rudnicki and Rice, 1975; Vardoulakis and Sulem, 1995). Crack propagation in rock features a complex behavior that usually needs to be studied numerically. The difficulty with the numerical rock softening and damage analysis is the dependency of the failure load on the mesh resolution. For example, it is well known that in the finite element analysis of rock softening, by reducing the mesh size, the dissipated energy in the fracture zone or shear band reduces which is not physically meaningful (Bazant and Planas, 1998). To remedy this difficulty, non- local continuum (Bazant and Pijaudier-Cabot, 1988), Cosserat continuum (De Borst and Muhlhaus, 1991) and gradient plasticity techniques (De Borst and Muhlhaus, 1992) have been proposed in the literature. Numerical continuum methods based on micromechanics and flaw distribution have been proposed as well. In a multi-laminate model (Pande and Xiong, 1982), dis- crete failure planes are considered in otherwise a continuum system. The original multi-laminate model later on was renamed a microplane model in that a geometric damage tensor was introduced (Carol et al., 1991). In a microplane model, a non-local averaging technique needs to be implemented to avoid the sensitivity of the results to mesh resolution. A displacement discontinuity technique was used by Van de Steen et al. (2003) to model rock fracture. This technique was successfully implemented in simulating rock spalling around an underground excavation, although no comment qu uniaxial compressive strength R cylinder radius Rave average radius of cylinders Rmax maximum radius of cylinders Rmin minimum radius of cylinders sb shear bond of a contact t specimen length in a Brazilian test �"axial axial strain increment �"lateral lateral strain increment ��axial axial stress increment � friction angle � friction coefficient of a contact �0 genesis pressure �1 axial stress and maximum principal stress �3 confining pressure and minimum principal stress �t Brazilian tensile strength � Poisson’s ratio �0 apparent Poisson’s ratio 194 A. Fakhimi and T. Villegas was presented regarding calibration of the model. Fractures can be studied as math- ematical cracks as well. In this case, linear fracture mechanics or cohesive crack models (Bazant and Planas, 1998) can be used, although extension of these models to biaxial compressive loading is not straightforward. In the discrete element approach of rock fracture, the interaction of particles is con- sidered explicitly in the model. This provides an embedded internal length in the for- mulation of the problem that prevents the numerical difficulties observed in the classical local continuum softening models. The application of this technique in geomechanics was initially proposed by Cundall and Strack (1979) and later on extended and implemented in the PFC computer programs (Itasca, 1999). With this technique, crack propagation, rock softening and dilation, anisotropic rock behavior in damaged zone, cyclic energy dissipation, spontaneous emission of acoustic energy and Coulomb failure envelope can be reproduced. Nevertheless, this powerful technique has a few shortcomings. Firstly, it needs calibration, i.e. some micro-mechanical properties must be specified to result in a synthetic material with desired macro-mechanical properties such as Young’s modulus, Poisson’s ratio, tensile strength, and failure envelope. For a simple contact bond model implemented in this paper, the micro-mechanical properties for interaction of two circular particles or cylinders are normal and shear spring stiffnesses (kn, ks), normal and shear bonds (nb, sb) to glue the particles in contact and Coulomb friction coefficient (�) which is activated when the bond between any two particles is broken. Therefore, these five parameters must be obtained during the calibration step. The calibration can be performed by using a trial and error approach in that extensive effort is needed to acquire appropriate micro-parameters to result in macroscopic properties similar to those of a specific rock. A more efficient method was employed by Huang (1999) for calibration of a circular particle interaction (CPI) model in that a dimensional approach was established and some dimensionless parameters were introduced. The second problem with a CPI model is the difficulty of calibrating it for the failure envelope. In general, a CPI model shows a friction angle which is smaller than that for a hard rock (Potyondy and Cundall, 2001). To increase the friction angle of a particle assemblage, non-circular particles can be used (Ting et al., 1993), but this requires additional computational time to detect contact of adjacent particles; contact detection of circular particles is much faster compared to that of other particle shapes. The third difficulty of a CPI model is ratio of unconfined compressive strength to tensile strength which is usually lower than that of a rock. Fakhimi (2004) proposed a slightly overlapped circular particle interaction (SOCPI) model to resolve the above second and third problems. The calibration issue of a SOCPI material that is based on the dimensional analysis is addressed in this paper. The calibration graphs introduced were used to develop a synthetic model for Pennsylvania Blue Sandstone. Several physical and numerical tests are reported to verify the validity of the numerical model. All the numerical analyses in this paper were conducted using CA2 computer program which was developed by the first author (Fakhimi, 1998). CA2 is a hybrid finite difference-discrete element code for large deformation static or dynamic de- formation of solid materials. The first step in discrete element modeling of rock, using CA2, is sample preparation. This step involves with generation of walls that define the domain of analysis. The walls in CA2 are not rigid but are discretized to smaller elements and are used to confine the cylinders. Next, the cylinders with uniformly Dimensional Analysis in Calibration of a Discrete Element Model 195 distributed radii are generated and inflated by increasing the cylinders radii. An explicit algorithm is used to obtain the equilibrium configuration of this particle assembly. At this stage of analysis, a contact bond model with no friction and bond at the contact points is assumed for the interaction of cylinders. After equilibrium, the induced stress in the surrounding walls is examined. By changing the number of cy- linders, the induced stress can be modified. This stress is called genesis pressure (�0) and its role is to create a small overlap between cylinders in contact. This small overlap is responsible for obtaining more realistic friction angle and ratio of compres- sive to tensile strength of the synthetic material (Fakhimi, 2004). The last step in sample preparation is to initialize the normal and shear contact stresses to zero and to introduce normal and shear bonds and friction between the cylinders in contact. More details about sample preparation are given in Fakhimi (2004). 2. Calibration Charts Calibration of a synthetic model made of rigid circular cylinders is a tedious task if it is undertaken by using a trial and error approach. A much more convenient and efficient method is to introduce some calibration curves using a dimensional analysis. In this paper, a synthetic material was intended to be calibrated for Young’s modulus, Poisson’s ratio, failure envelope and tensile strength. A SOCPI model is calibrated by assigning appropriate parameters including normal and shear spring constants (kn, ks), normal and shear bonds (nb, sb), friction coefficient (�), genesis pressure (�0), average cylinders radius (R) and specimen width (D). The parameters that dictate the elastic properties of a SOCPI material are kn, ks, �0, R and D. These five parameters have two independent dimensions, i.e. force (F) and length (L). Hence, considering the p theorem in dimensional analysis (Sedov, 1993), 5� 2¼ 3 dimensionless parameters are involved in describing the elastic deformation of the synthetic material, i.e. E0 ¼ kn f ðR=D; ks=kn; �0=knÞ; ð1aÞ �0 ¼ gðR=D; ks=kn; �0=knÞ; ð1bÞ where f and g are dimensionless functions and E0 and �0 are apparent Young’s modulus and Poisson’s ratio of the synthetic material, defined in a uniaxial test as: E0 ¼ ��axial=�"axial ð2aÞ �0 ¼ ��"lateral=�"axial: ð2bÞ If a numerical experiment is considered as a biaxial plane strain test, true Young’s modulus and Poisson’s ratio of synthetic material can be obtained from the following equations: � ¼ �0=ð1 þ �0Þ; ð3aÞ E ¼ E0ð1 � �2Þ: ð3bÞ Several numerical uniaxial tests were conducted on the synthetic material with R=D¼ 0.013 and different ks=kn and �0=kn values, using CA2 computer program. 196 A. Fakhimi and T. Villegas The numerical uniaxial test set up together with the synthetic specimen is shown in Fig. 1. The upper and lower platens or walls are made of finite difference grid and are frictionless. The upper wall is fixed vertically while the lower one was moved with a constant velocity to compress the specimen. Two small and very flexible walls were glued to the lateral sides of the specimen. These walls were like measurement devices to ac- Fig. 1. Numerical unconfined compressive test set-up Fig. 2. Elastic constants of synthetic material versus ks=kn for different �0=kn, a normalized apparent Young’s modulus, b apparent Poisson’s ratio Dimensional Analysis in Calibration of a Discrete Element Model 197 quire average lateral deformation of the middle of the specimen and consequently the Poisson’s ratio. The results of the numerical analyses are shown in Fig. 2a and b. Para- metric studies and also the work by Potyondy and Cundall (2001) show that E0 and �0 are almost independent of R=D values provided that this ratio is kept small enough. Therefore, by introducing the macroscopic parameters E0 and �0 and assuming a value for �0=kn, normal and shear spring constants can be estimated from Fig. 2a and b. The uniaxial strength of the synthetic material (qu) was obtained by conducting several numerical tests. This parameter is a function of �0, kn, ks, R, D, nb, sb, and �. These eight parameters have two dimensions. Hence, dimensionless unconfined com- pressive strength can be expressed as: quR=nb ¼ f1ð�0=kn; ks=kn; knR=nb; nb=sb;R=D; �Þ; ð4Þ where f1 is a dimensionless function. Numerical analyses demonstrate that qu is a weak function of � and kn R=nb (Villegas, 2004). In addition, it has been shown that the biaxial compressive strength of a synthetic material made of circular cylinders is not a function of R=D value, provided that R=D is kept a small number (Potyondy and Fig. 3. Dimensionless uniaxial strength versus ks=kn for a �0=kn¼ 0.007 and b �0=kn¼ 0.08 198 A. Fakhimi and T. Villegas Cundall, 2001). Thus, only three dimensionless parameters namely, �0=kn, ks=kn and nb=sb were needed to consider. By conducting uniaxial tests on the synthetic materials with R=D¼ 0.013 (Rmin¼ 0.44 mm and Rmax¼ 0.58 mm with a uniform random dis- tribution of cylinders radii), �¼ 0.5, kn¼ 50 GPa, and nb¼ 20,000 N=m, two charts, each for a specific value of �0=kn (0.007 and 0.08) were obtained. Figure 3a and b show these charts where the dimensionless uniaxial strength is plotted versus ks=kn for different nb=sb values. Similar charts for �0=kn¼ 0.04 and 0.1 are reported in Villegas (2004). Comparison of these figures verifies that the dimensionless unconfined com- pressive strength of synthetic material increases as higher values of genesis pressure is used. This should be expected, as with a higher genesis pressure, the number of con- tact points per unit area of the specimen increases which results in a stronger speci- men. It is interesting to observe the effect of ks=kn parameter on the specimen strength. Figure 3 indicate that by reducing the ks=kn parameter to a value less than 0.5, the dimensionless strength decreases. The crack pattern, in the synthetic material is also affected by the ks=kn ratio. This is demonstrated by the crack pattern for the synthetic specimens with different ks=kn value (Fig. 4). A crack in CA2 is a line perpendicular to the line connecting centers of cylinders with broken bond. It is evident from Fig. 4 that for a small ks=kn value, a diffuse crack pattern is obtained, in contrast to a shear banding when a greater value of ks=kn is used. The reason for this observation is that by increasing ks value, greater microscopic shear stresses and hence shear cracks are developed in the specimen with the consequence of shear banding. In addition, it can be shown that a synthetic specimen with a small ks=kn value, features a substantial amount of hardening in its stress-strain diagram, while a greater value of ks=kn, results in a more brittle behavior (Villegas, 2004). Several numerical biaxial tests were conducted to evaluate failure envelope of synthetic material. Strength of synthetic material (�1) is a function of �0, R, D, �, ks, Fig. 4. Crack pattern in the damaged specimen with a ks=kn¼ 0.05, b ks=kn¼ 0.55, and c ks=kn¼ 1.5 Dimensional Analysis in Calibration of a Discrete Element Model 199 kn, nb, sb, and the confining pressure (�3). Similar to uniaxial strength, �1 is a weak function of � and knR=nb. Therefore, the biaxial dimensionless strength can be written as: �1R=nb ¼ g1ð�3R=nb; ks=kn; �0=kn; nb=sbÞ; ð5Þ where g1 is a dimensionless function. The normalized failure envelopes for dif- ferent values of �0=kn are shown in Fig. 5a–c. In this series of numerical anal- Fig. 5. Normalized failure envelopes for a �0=kn¼ 0.007, b �0=kn¼ 0.04, and c �0=kn¼ 0.08 200 A. Fakhimi and T. Villegas yses, the following parameters were fixed: R=D¼ 0.013 (with Rmin¼ 0.44 mm and Rmax¼ 0.58 mm), �¼ 0.5, kn¼ 50 GPa, ks¼ 27.5 GPa and nb¼ 20,000 N=m. Different values of nb=sb (0.05, 0.1, 0.2, 0.3, 0.4, 0.5) were used in the analyses. Figure 5 indicate that by increasing the genesis pressure, specimens with higher strength are obtained. In addition, it is noticed that with larger �0=kn values together with smaller values of nb=sb parameter, curvature of failure envelope similar to that of a rock can be simulated. The slope of failure envelope in a �1, �3 space is equal to tan 2(45þ�=2) in that � is the friction angle. Inspection of the curves in Fig. 5 for higher �0=kn and smaller nb=sb values show that friction angles of 60 � or more can be simulated. This has become feasible due to the application of non-zero genesis pressure. To investigate the role of microscopic friction coefficient (�) on the failure envelope, two sets of biaxial tests with �0=kn¼ 0.08, ks=kn¼ 0.55, nb=sb¼ 0.2, kn R=nb¼ 1269, and R=D¼ 0.013 were conducted. For the first set of biaxial tests, a friction coefficient of 0.5 was used while for the second set, � was 0.8. Figure 6 shows the peak and residual failure envelopes for these two synthetic materials. This figure demonstrates that while � has a small effect on the peak failure envelope, its role on the residual strength is negligible. Consequently, a great percentage of the induced friction, in the numerical model, is due to the shear displacement along the developed irregular failure surface and to a much less extent the microscopic friction angle (�) plays a role. Some numerical Brazilian tests were performed to evaluate tensile strength of synthetic material. The numerical test set up for the Brazilian test is shown in Fig. 6. Normalized peak and residual failure envelopes for two different friction coefficients Fig. 7. Numerical Brazilian test set-up Dimensional Analysis in Calibration of a Discrete Element Model 201 Fig. 7. The upper platen is fixed in the vertical direction while the lower platen was moved with a constant velocity of 0.2� 10�8 m=step in the upward direction. Both upper and lower platens were discretized to the smaller finite elements. Two series of numerical Brazilian tests were conducted. In the first series, for specific values of �0=kn, the tensile strength was obtained by changing ks=kn and nb=sb parameters. The tensile strength was obtained from the following equation (Goodman, 1989): �t ¼ 2p=�Dt; ð6Þ where p is the maximum applied load, D is the specimen diameter and t is its length that is equal to unity in the numerical tests. Figure 8a and b show the dimensionless tensile strength (�t R=nb) versus ks=kn parameter for �0=kn¼ 0.1and 0.04, respectively. Similar to the observation with the unconfined compressive strength, it is evident that increase in genesis pressure results in higher tensile strength. The reason for the observed maximum points of dimensionless tensile strength as shown in Fig. 8a is not clear to the authors. One possible explanation is that by increasing nb=sb parameter with higher values of �0=kn, the intensity of micro shear cracks are increased in the model that results in change of the trend of the failure load. Fig. 8. Dimensionless Brazilian tensile strength versus ks=kn for a �0=kn¼ 0.1 and b �0=kn¼ 0.04 202 A. Fakhimi and T. Villegas In the second series of Brazilian tests, the size effect was investigated. Figure 9