Indian Geotech J https://doi.org/10.1007/s40098-018-0304-8
ORIGINAL PAPER
A Study on Crack Initiation and Propagation in Rock with Pre-existing Flaw Under Uniaxial Compression G. Sivakumar1 • V. B. Maji1
Received: 16 November 2017 / Accepted: 6 April 2018 Ó Indian Geotechnical Society 2018
Abstract A study for crack initiation and propagation in rock with a pre-existing flaw is carried under uniaxial compression. The aim is to understand the influence of crack initiation stress and peak stress for narrow flaws in rock by conducting laboratory experiments and subsequent numerical simulations. Specimens are prepared using Gypsum material with the incorporation of flaws at varying flaw angles from 15° to 75° having three different strengths. The results of the crack initiation and peak stresses for different flaw angles along with crack growth patterns are observed. Numerical simulations using ABAQUS adopting extended finite element and cohesive zone model is attempted to capture the crack initiation and propagation in rock specimens. For closed crack condition, the numerical model considers frictional behaviour across the flaw surface and model could capture the crack initiation and peak stresses along with propagated crack patterns irrespective of the material strengths. Both laboratory and numerical studies reveal that the crack initiation stress and peak stress gets initiated sooner for the flaw having 45° when compared to all other angles which are in contrast to open flaw conditions. Also at this angle, a stable curvilinear wing cracks are observed irrespective of the material strengths. In laboratory experiments, pre-dominant wing cracks are observed for a material having high strengths that influence the final pattern in specimens while in numerical analysis, all the specimens shows clear visible wing crack propagation and leads to failure of the material.
& G. Sivakumar
[email protected] 1
Department of Civil Engineering, Indian Institute of Technology Madras, Chennai 600036, India
Keywords Narrow flaws Extended finite element (XFEM) Cohesive zone model (CZM) Crack initiation Peak stress
Introduction Rocks inherently exhibits anisotropy due to the existence of discontinuities such as joints, fractures, faults, and bedding planes. The existence of these discontinuities can alter the state of stress in the rock mass when compared with corresponding intact rocks. When such rockmass containing fractures are subjected to loading, new cracks may develop altering their strength and deformation characteristics. These newly formed cracks will propagate towards the direction of the loading and may lead to macroscopic failures. The fracture mechanics theory can convincingly describe the effect of those crack formation and corresponding growth under various loading. Over the past few decades, many researchers attempted to study the rock fracture behaviour based on crack initiation stress and their propagation [1–8]. The Fracture mechanics concept was first introduced by Griffith [1] proposing a fracture initiation criteria for brittle materials by considering a single elliptical open flaw. Griffith developed the criteria based on the local stress state for an elliptical open hole in an infinite plate derived by Inglis [2]. The Griffith theory assumes that the existing fracture initiates once the stress at the crack tip exceeds the molecular cohesive strength of material [3] or the uniaxial tensile strength of the material [4]. However there are some limitations, Hoek and Bieniawski [5] found that Griffith’s fracture criteria suitable for crack initiation but does not provides the information on the direction of crack propagation. Also, Griffith’s criteria
123
Indian Geotech J
are developed only for pure brittle materials where plasticity is not considered. But in reality, a fracture process zone exists in front of the pre-existing flaw tip where the plastic deformation occurs [6]. Based on analytical models by Bareblatt [6] and Dugdale [7], Hillerborg [8] proposed a new model known as cohesive zone model (CZM) which considers the plasticity in linear elastic fracture mechanics (LEFM). Some researchers have conducted experimental studies on crack growth in rock or rock like materials subjected to compressive loads [9–15]. Bobet and Einstein [9, 10] and Sagong and Bobet [11] have done an investigation on gypsum, a rock like material for capturing crack pattern and their propagation for open and narrow flaws. A better visualization of crack propagation types was done by Wong and Einstein [12, 13] in their study by using the high speed camera. Though the crack initiation stress is covered in these previous studies, the effect of friction behaviour was not considered in their experimental study [14, 15]. Also the fracturing process was found to be similar for both open and closed flaws, the effect of stresses found to be varying. With advancement in computational techniques, the studies on crack stresses and its propagation are extended into numerical methods [8, 16–19]. There are various numerical methods like finite element method (FEM) [16–19], discrete element method (DEM) [20], hybrid finite-discrete element methods (FDEM) [21, 22] and displacement discontinuity method (DDM) [23] to study the crack initiation and propagation by different criteria based on fracture mechanics. FEM is considered to be one of the most widely adopted methods for modelling the material having a pre-existing flaw in a continuum medium. In 1976, Hillerborg [8] made an attempt to consider CZM into FEM for concrete materials. Later it was adopted by Reyes and Einstein [16], and Gonc¸laves da Silva and Einstein [18] on FEM for rock and rock like materials by considering three different fracture criterion namely principal stress, principal strain and energy criterion. Though the obtained results agree reasonably well with experimental studies, the numerical model does not account for frictional effect between the flaws. Most of these studies considered rock with open flaws, subjected to compressive loading where the flaw remains open throughout the loading history. Rock domain may exhibit both open flaws and narrow (closed) flaws as shown in Fig. 1, where the latter gets closed during loading leading to surface friction between the flaws [19]. Due to this behaviour there can be changes in crack initiation and peak stresses of the material. The studies on the crack having narrow flaw subjected to compressive loading revealed that the crack initiation stress was found to be higher than in open flaw [14–16]. The influence of frictional behaviour across the flaw surface is accounted into analytical studies by Steif [24], Horii and Nemat-
123
Nasser [25], Ashby and Hallam [26] and Baud et al. [27]. Xie et al. [19] carried out the numerical analyses to understand the effect of crack surface friction using extended finite element method (XFEM) based on LEFM. However, their adopted numerical criteria considered only the crack initiation stress. The post influence of propagation after crack initiation and the material threshold strength (peak stress) are not covered. And it is important to have an appropriate criterion that can determine crack initiation stress with respect to peak stress for engineering design and practice. The present study is an attempt to determine the peak stress after crack initiation in the rock like material having narrow flaw under uniaxial compression. An experimental study is carried out in gypsum having a narrow flaw with three different material strengths. The material is subjected to uniaxial compressive loading, where the crack growth behaviour till the peak strength is studied. The flaw angles are varied from 15° to 75° to understand the phenomenon of the crack initiation stress and the peak stress of the material. Obtained results are subsequently validated with a numerical model. The study uses commercially available FE tool ‘‘ABAQUS’’ to model using CZM with XFEM. The major advantage of XFEM is it does not require mesh refinement in front of the existing flaw. CZM model adopts traction separation law, a damage criterion used for describing the crack initiation and propagation process. In the present study, both experimental and numerical studies are conducted and their results in terms of both crack initiation stress and peak stress are discussed. Experimental Studies on Crack Growth Several researchers conducted experimental studies to understand the behaviour of crack propagation subjected to uniaxial compression for different materials such as Hoek and Bieniawski [7] on glass, Lajtai et al. [28] on Plaster of Paris, Bobet and Einstein [10], Sagong and Bobet [11], Park and Bobet [14, 15] on gypsum, Wong and Einstein [12, 13] on gypsum and marble. Most of the previous studies considered the flaw to be both open and closed flaws where the study was limited to understand the crack growth and coalescence pattern between the multiple flaws. The staged behavior of a single flaw in closed condition subjected to uniaxial compression was first explained by Lajtai [28]. Laboratory tests on Plaster of Paris material, having specimen dimensions 300 9 600 9 600 (76 mm 9 152 mm 9 152 mm) with flaw aperture width 0.0200 (0.05 mm) of variable lengths 0.500 , 100 , 1.500 and 200 was carried out. Lajtai found that before the crack gets initiated the material will behave in an elastic manner. Later, new cracks are formed which can be of tensile and shear natures. After initiation, the material yields its threshold value
Indian Geotech J
Fig. 1 Open and closed flaw as commonly found in a rock mass
followed by formation of the shear zone around the crack tip. Studies with various rock-like materials commonly observed two types of cracks namely primary and secondary crack (Fig. 2). Primary wing cracks are in tensile nature or tensile shear nature and secondary cracks are shear in nature which appears after the initiation of primary cracks. A primary wing crack starts at the tip of the flaw and propagates in a curvilinear path under loading. Secondary cracks appear in a quasi-coplanar and oblique direction around the crack tip along with shear zones. The secondary cracks not always being observed [19], unlike primary cracks which were found to be dominating. With advanced experimental techniques using a high-speed camera, Wong and Einstein [13] observed four different types of primary cracks under uniaxial compression as shown in Fig. 3.
Numerical Study To capture the crack growth according to fracture mechanics, the present study uses cohesive zone model (CZM). The model is based on the concept of damage evolution around discontinuity (pre-existing flaw) that occurs after crack initiates. CZM uses traction separation law for modelling the discontinuity. In the present study, the adopted model is combined with XFEM approach using FE tool ‘‘ABAQUS’’. The basic concepts of XFEM and CZM modelling are explained in this section. Cohesive Zone Model The CZM follows a traction–separation relation across the elements at flaw tip to capture the damage around the discontinuity. The behaviour of the material is assumed linear elastic prior to the initiation of the crack from existing flaw [29]. When the load is applied to the material, the model assumes that there exist the traction components in the existing flaw. The traction components present failure in tension (Mode I) or by shear in-plane (Mode II) or by out-plane shear (Mode III) or a mix of all three modes. For a two dimensional problem, nature of the traction components in tension (Mode I) and shear in-plane (Mode II) failure modes are generally considered (Fig. 4). Crack Initiation and Propagation There are various crack initiation criteria available for rock or rock like materials namely maximum principal stress criteria, maximum principal strain criteria and maximum shear stress criteria. The present study adopts widely used maximum principal stress criterion (f) for describing the onset of the crack initiation,
Fig. 2 A simplified crack pattern in pre-existing flaw subjected to uniaxial compression
123
Indian Geotech J Fig. 3 Crack propagation observed from a narrow single flaw of molded Gypsum [13]. a Type 1 tensile crack (tensile wing crack), b Type 2 tensile crack, c Type 3 tensile crack and d mixed tensile-shear crack
an oblique direction with respect to the crack plane (Fig. 5).
Fig. 4 Traction components acting on pre-existing flaw subjected to compression
Damage Evolution After Crack Initiation As per traction–separation law when the load is applied, traction components acting at the flaw will gradually increase towards the maximum material limit in tension, Tn and shear, Ts (Fig. 6). Once the load exceeds the material limit, a new crack will get initiated with decrease in traction components at the flaw tip [30],
f ¼
hrmax i Tn;s
ð1Þ
The relation denotes that the discontinuity gets initiated when traction components Tn,s present across exceeds hrmaxi, the assigned maximum principal stress value. The symbol hi represents the Macaulay bracket, signifying that the damage does not initiate for a purely compressive stress state. When the maximum principal stress ratio (Eq. 1) approaches a value of unity the damage initiates. As per criteria, the new initiated crack under mixed mode (considering both tension and shear traction) will propagate in
Fig. 5 Propagated new crack based on mixed mode condition
123
tn ¼ ð1 DÞ Tn ; where Tn ¼ Kp dn
ð2Þ
ts ¼ ð1 DÞ Ts ; where Ts ¼ Kp ds
ð3Þ
where tn and ts are instantaneous tractions in tension and shear respectively, Tn and Ts are a maximum limit of tractions respectively. dn and ds are the relative displacements due to tension and shear tractions respectively acting at the pre-existing flaw. The relative displacements also termed as the displacement
Fig. 6 Traction–Separation law
Indian Geotech J
discontinuity which measures the discontinuity (preexisting flaw) deformation due to tension (Mode I) and shear (Mode II) tractions with reference to the discontinuity plane. The term Kp is the penalty stiffness which controls the amount of penetration between the elements in the contact region of flaw surface when subjected to closure [21, 22]. If the contact between the flaw surface is assumed to be open, then their value becomes infinity. If smaller stiffness value is adopted, it will lead to the penetration of the elements. When the cracks are initiated, the tractions get degraded from its threshold limit represented by damage variable, D. It will have a value of zero when a crack initiates and it monotonically becomes unity upon further loading leading to complete degradation. A linear damage model is used to describe the damage variable which is given by: 0 dfm dmax m dm D ¼ max f ð4Þ dm dm d0m where dfm is the mixed mode relative displacement that is observed at the complete failure. dmax is the maximum m value of the mixed mode relative displacement attained at that loading history and d0m is the mixed mode relative displacement that is observed at the crack initiation. When dmax reaches d0m , a new crack gets initiated where D m is zero. The initiated crack will be propagating, when f d0m \ dmax m \ dm and D increases linearly from the value of f zero to one. Once dmax m [ dm the material attained its full damage where D reaches the value of one. In the present study, relative displacement is used in conjunction with mixed mode function to represent the degradation of traction components under both tension and shear along with mode mix function. Relative displacement under mixed mode condition (dm), also termed as discontinuity displacements was first introduced by Camanho and Davila [30]. It correlates separation of the discontinuity due to tension and shear stress components acting on it from the threshold limit of traction components until its complete degradation and is expressed as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dm ¼ hdn i2 þd2s ð5Þ
u¼
2 Ts tan1 p Tn
ð6Þ
The propagated crack grows along the element boundary which is mainly dependent on the mesh geometry assigned around the discontinuity of the material. Using XFEM, the crack path grows arbitrarily into the elements without mesh refinement. Interface Friction When the material having a very narrow gap flaw is subjected to compression, they experience friction at their interface. In the present study, the pre-existing flaw at both upper and lower inner surface comes in contact, where the interaction between those surfaces contributes friction. The current model considers a friction coefficient parameter for sliding movement between the flaw surfaces after its contact. Usually, the sliding friction coefficient has a lower value than the static friction (to resist the movement of the slide). In the model, the flaw gets to slide at a friction stress (ss) which is in terms of sliding friction coefficient (f) and normal and shear traction. It is given by the relation as follows, ss ¼ ts tn f
ð7Þ
Basic Concept of XFEM Modelling The extended finite element method (XFEM) was first introduced by Belytschko and Black [31] for applications exhibiting crack propagation and is an extension of FEM. This method allows a new crack to get propagated without remeshing at the discontinuity (flaw) tip. Thus making it suitable for fracture process analysis where mesh does not require a structure to initiate a crack and further crack path definitions. In XFEM, the elements around the fracture surface are enriched with functions having additional degrees of freedom based on the partition of unity (Fig. 7). The enrichment functions represent the jump in displacement across the discontinuity surfaces and its
The mode mix function is used in conjunction with relative displacement in the numerical analysis. The mode mix function consists of tension and shear traction components. The newly developed crack will propagate in the oblique direction with respect to existing discontinuity plane and further the crack propagated towards the loading direction. The mode mix function is expressed as: Fig. 7 Discontinuity present across the element domain with Heaviside enrichment
123
Indian Geotech J
approximation for a displacement vector function with the partition of unity enrichment is u¼
N X
NI ðxÞ½uI þ HðxÞaI
ð8Þ
I¼1
where NI(x) is nodal shape functions, uI is nodal displacement present in all the nodes and I = N represents the number of gauss integration points. aI is the enriched nodal degrees of freedom vector and H (x) is called Heaviside enrichment across the crack surface (Fig. 4) which is given by, 1 if ðx xÞ n 0 HðxÞ ¼ ð9Þ 1 otherwise where, x is an integration point, x* is the point on the crack closest to x, and n is the unit normal to the crack at x* (Fig. 8). In the present numerical analysis, nodes related to the elements that have discontinuity are enriched using XFEM. Fig. 9 a Mold for sample preparation and b cast gypsum specimen with a single flaw orientated by 45° (b)
Methodology Experiment Study Sample Preparation The experimental studies are carried out with prepared specimens with dimension 76 mm (width) 9 152 mm (height) 9 30 mm (thick) of gypsum plaster with three different compressive strengths of 15, 35 and 60 MPa representing low, medium and high strengths respectively. A flaw of length 12.5 mm made at different inclinations similar to some of the earlier experimental studies [10, 13]. Commercial grade gypsum powder is mixed with a predetermined volume of water and the resulting slurry is poured into the fabricated mold (Fig. 9a) to prepare the specimens. For creating a narrow flaw, a metal shim of a very smooth surface having a thickness of 0.1 mm is used, in mid-section of the specimen mold. It is introduced before pouring the slurry as shown in Fig. 9. The specimen is allowed for curing and after some time, the metal shim is removed before gypsum attains its full strength. Specimens are prepared with five different flaw angles varying from 15° to 75° (with 15° interval) with three different material
strengths. The specimens are kept for complete curing for 24 h by air drying at room temperature and subsequently kept in the oven at a constant temperature of 40 °C for 4 days [10, 13, 15]. The physical and mechanical properties of the material obtained after 5 days of preparation are listed in Table 1. Prefix L/M/H are used to denote the three different material strength of each specimen respectively with the flaw angle inclination in degrees. Sample Testing All the prepared gypsum specimens of three different compressive strengths are tested in compression testing machine having a loading capacity of 500 kN (Fig. 10). The load is applied from bottom platen along the height of the specimen at a slow rate of 0.031 mm/min while the upper platen is fixed [10, 13]. Crack initiation load and the respective peak load are monitored from the load cell fixed at upper platen. The final failure pattern after it attains peak strength is also captured for all the specimens. Numerical Modelling Determination of Input Parameters
Fig. 8 Illustration of a gauss point for jump function used across the crack surface
123
The input parameters are obtained using laboratory tests namely uniaxial compression, Young’s modulus and Poisson’s ratio as suggested by ISRM [32, 33]. Tensile strength is determined using Brazilian tensile test [34]. For CZM parameters, the samples are prepared and tested based on Centre Straight Notch Brazilian Disc (CSNBD)
Indian Geotech J Table 1 Input material parameters from laboratory experiments used for numerical analysis Material parameters
Low strength
Medium strength
High strength
Young’s modulus (GPa)
0.898
1.927
3.985
Poisson’s ratio
0.39
0.34
0.32
Tensile strength (MPa)
1.11
2.14
3.61
Tensile stress for mode I (MPa)
1.087
1.880
3.294
Shear stress due to mode II (MPa)
1.560
2.515
4.946
Relative displacement under mixed mode (lm)
8
0.72
0.45
Penalty stiffness (GPa/m)
4.49
9.64
19.95
Fig. 10 Uniaxial compression testing (CT) machine where prepared specimens are tested
fracture tests. Based on the procedure by ISRM [36], Atkinson et al. [37] and Al-Shayea [38], the tension and shear tractions across the flaw surface and their corresponding relative displacements for both Mode I (pure tension) and Mode II (pure shear) failures are found. The relative displacement parameter under mixed mode condition is calculated based on Eq. 4. It is observed from the Table 1 that the relative displacement parameter exhibits relatively higher value for low strength indicating more plasticity at the flaw tip. Based on laboratory test procedure [35] for gypsum having smooth and planar surface the sliding friction coefficient value found to be varying from 0.23 to 0.28 with an average value of 0.25. For normal behaviour, a stiffness parameter is used for controlling the penetration of flaw surface when subjected to loading. The contact penalty stiffness is chosen five times that of Young’s modulus of the material based on the recommendation by Mahabadi [21] and Lisjak et al. [22].
integration method and the corresponding mesh is shown in Fig. 11. The boundary conditions of the model are made in
Numerical Model in ABAQUS The element type used for this model is CPS4R, a 4-node bilinear plane stress quadrilateral and solved by reduced
Fig. 11 A numerical model of gypsum specimen in ABAQUS with a flaw at 45° angle
123
Indian Geotech J
L15
L30
L45
L60
L75
M15
M30
M45
M60
M75
H15
H30
H45
H60
H75
Fig. 12 Final failure pattern of fully developed crack growth observed for three different material strengths with Type 1 (T1), Type 2 (T2) and Type 3 (T3) failure modes based on Wong and Einstein [13]
accordance with the laboratory tests, where the top end is restrained in the Y-direction and the bottom end is
123
subjected to loading. Movement in the X-direction is not restricted so as to allow free crack propagation. The
Indian Geotech J
analysis is carried out in dynamic implicit quasi-static procedure and load is applied in a displacement controlled manner with a loading rate of 0.031 mm/min as applied during laboratory studies.
Results Experiment Results The final failure patterns for all tested specimens are shown in Fig. 12. Specimens subjected to uniaxial loading exhibit predominant wing failure at both the tips of the flaw during the initial stages of the crack development. The primary wing cracks are observed to initiate in the oblique plane direction of the flaw surface. Cracks gradually orient towards the direction of loading following a curvilinear path. The wing cracks are found to be more prominent for the higher strengths. Specimens with flaw angle varying from 15° to 60° are found to have some spalling at the crack tip region. This is observed commonly irrespective of the material strength. But for specimens having a flaw angle of 75°, spalling is observed away from the crack tip. Loading of the specimen leads to closure of flaw and subsequent spalling induces high concentration of compressive stress at the flaw tip. When the flaw angle increases towards the loading direction, the frictional resistivity between the flaws reduces. Therefore the specimens with steeper flaw angles (flaw angle [ 60°) are less prone to spalling at the flaw tips. Similar observations are reported by Xie et al. [19] in their numerical study. No subsidiary cracks are formed on the specimen apart from wing cracks. But when the material strength increases additionally subsidiary cracks are observed specifically in high strength specimens. It can be assumed when the material strength increases, the shear resistance between the propagated crack surfaces causes more subsidiary cracks along the path of main crack propagation [39]. Further application of load results in complete rupture of the specimen. The crack initiation stress and the peak stress for different flaw angles (b°) are shown in Fig. 13. Flaw angle of 45° is found to have less crack initiation stress as well as peak stress for all the material strengths. The variation of crack initiation stress and peak stress with respect to different flaw angles is found to be similar. The difference between the crack initiation stress and peak stress is observed to be the least for high material strength. The peak stress occurs immediately after crack initiation stress for the flaw angles 15° and 75°, irrespective of the material strength. Figure 14 shows the plots of the stress and strain observed in the vertical direction for all flaw angles in medium strength material. The crack initiation stress is
Fig. 13 Crack initiation and peak stress observed in experimental studies
Fig. 14 Laboratory stress versus strain plot observed for medium strength gypsum specimens with angles from 15° to 75°
identified as the stress corresponding to the sudden dip observed in the stress–strain plot. Specimens attain failure in small strain for the flaw angle of 45°, where the wing cracks are found to have a longer curvilinear path when compared to other angles. The cracks are initiated at larger strains for the flaw angles of 15° and 75°. Numerical Results The failure pattern observed in the numerical analyses [40] with flaw angles varying from 15° to 75° is shown in Fig. 15. The wing crack growth is clearly seen to propagate
123
Indian Geotech J
L15
L30
L45
L60
L75
M15
M30
M45
M60
M75
H15
H30
H45
H60
H75
Fig. 15 Failure pattern observed from ABAQUS for all three strengths with flaw angles varied from 15° to 75° where Type 1 (T1) and Type 2 (T2) failure modes of Wong and Einstein [13]
123
Indian Geotech J
from the tips of the flaw to the material boundary in the numerical analysis. Additional cracks are formed especially in high strength material after the model reaches its threshold limit. The specimens do not have significant subsidiary cracks in models with a low strength. More subsidiary cracks are observed for medium and high strength materials similar to experimental results. The numerical model used in this study can determine both the crack initiation stress and the peak stress for all strengths (Fig. 16). The variation of crack initiation stress and peak stress with respect to different flaw angles is found to be similar and is in good agreement with the results of the experimental studies. Peak stress is reached immediately after the crack initiation stress, except for the model with 45° flaw angle. The curvilinear path is observed to be much more pronounced for the model having flaw angle of 45°. However, for the model of high strength material with flaw angle of 45°, the wing crack observed is relatively smaller compared to low and medium strength material. Moreover the difference between crack initiation stress and peak stress becomes smaller when the material strength increases. Based on the relative displacement i.e. input parameter for materials (refer Table 1), the higher strength material which has lower value of relative displacement implies less plasticity at the flaw tip. This causes the model to reach the material peak strength immediately after its crack initiation stress. Curvilinear wing path is relatively small for higher strength specimens as compared to low and medium strengths. This is due to the fact that as material strength increases the slippage between the flaw decreases. Therefore the
propagated crack will nucleate in unstable manner with formation of more subsidiary cracks along the crack path. Figure 17 shows the stress versus strain plot for medium strength model to understand the characterisation of stress with different flaw angles. A M45 specimen model is shown in Fig. 18a to demonstrate the propagation of crack based on CZM. Three points (P1, P2 and P3) are chosen along the wing crack path, where P1 is the point where crack initiates, P2 is the point where the curvilinear wing crack path orients towards the loading direction and P3 is the point where the propagated crack meets the boundary of the material. The maximum principal stress and the axial strain are monitored at the three points until the material fails (Fig. 18b). It is observed that at the point P1, crack gets initiated just before the axial strain of 1% and the maximum principal stress reaches a value of about 3.68 from 2.88 MPa. The maximum principal stress level drops down to zero at the axial strain value of 1% which is not clearly visible in Fig. 18b because the drop in the maximum principal stress occurs within a negligible change in the axial strains. This drop in the maximum principal stress indicates the development of wing crack propagation from the flaw tip. When the wing crack path reaches the point P2 at 1.03% axial strain, the maximum principal stress obtained is slightly lower than the previous maximum principal stress at point P1. The crack grows curvilinearly in a stable manner between these two strain levels. The wing crack path moves from point P2 towards the point P3 and attains a strain value of 1.187% where the model reaches its maximum threshold value. Formation of subsidiary cracks observed around the points P1, P2 and P3 leads to the
Fig. 16 Crack initiation and peak stress observed from numerical model
Fig. 17 Numerical stress versus strain plot obtained for medium strength with flaw angles varied from 15° to 75°
123
Indian Geotech J Fig. 18 Numerical model of M45 a points monitored for maximum stress evolutions and b plot showing maximum principal stress evolution observed with respect to strain
increase in the maximum principal stress value. The maximum principal stress at point P3 is found to have a higher value than that at points P1 and P2 where more number of subsidiary cracks are formed at the boundary in unstable manner. Further, loading leads to the complete failure of the material at an axial strain of 1.36% and the maximum prinicipal stress value drops after the formation of the subsidiary cracks. The numerical analysis sucessfully predict the crack propagation based on the adopted CZM where the maximum principal stress criterion have been utilised. Comparison Between Experimental and Numerical Studies It is observed from experimental and numerical studies that the initial cracks appear and propagate from flaw tip region exhibiting mainly Type 1, Type 2 and Type 3 tensile failure modes as reported by Wong and Einstein [13]. All the propagated cracks in the experiments are initiated in the oblique direction of flaw surface. Type 1 as well as Type 3 wing cracks are observed for H15 and H30 specimens in
123
the experiments but only Type 1 wing cracks are observed in the numerical analysis. Most of the specimens are found to have a combination of Type 1 and Type 2 cracks in the laboratory experiments. In the case of M30, M60, H30, H45 and H60 specimens, wing crack types are not fully observed due to the formation of spalling around flaw tip. Subsidiary cracks are also formed around the flaw tip in the numerical model along with the wing crack. Predominant Type 2 wing cracks are observed for flaw angle of 75° in both experimental and numerical studies. More subsidiary cracks are observed in both the experimental and the numerical studies, for high strength material with higher brittleness compared to materials of low and medium strengths. Figures 19 and 20 show the comparison of crack initiation stress and peak stress for both the experimental study and numerical analysis. Experimental values are found to be lower by 2–18% compared to the values observed in numerical analysis. This difference could be because numerical analysis is performed using two-dimensional models which may not be able to exactly simulate the three-dimensional nature of the physical experiments [20].
Indian Geotech J
two types are observed causing the wing crack to reach the material limit faster than experiments. However, the percentage of error between experimental and numerical results for the crack initiation stress as well as the peak stress for the high strength material is found to be comparatively lesser. It is observed from both experimental and numerical studies that the crack is found to initiate relatively earlier for the flaw angle of 45° compared to other angles, irrespective of the material strength. The crack initiation stress and peak stress therefore to have minimum value for the material having narrow flaw condition which is also reported by Ashby and Hallam [26].
Conclusions
Fig. 19 Comparison between the experimental and the numerical results of crack initiation stress observed for three different strengths specimens with respect to flaw angle
The crack initiation stress and corresponding peak stress in artificial rock specimens with a single narrow flaw under uniaxial loading are investigated with laboratory experiments and numerical analysis. The numerical model uses cohesion zone model (CZM) and is combined with XFEM in ABAQUS where the mesh does not require refinement around the discontinuity (flaw). This XFEM based numerical model could capture the cracks initiation and peak stress with primary wing crack and subsidiary crack formations. The numerical results are compared with the laboratory experiments and systematically verified. It was found that numerical values of crack initiation and peak stress and the crack growth pattern closely compare with the laboratory experiments. In both cases, the crack initiation and peak stresses do not increase with an increase in flaw angle, rather a minimum crack initiation stress and peak stress was found for a specimen having flaw angle of 45°. The differences between the crack initiation and peak stress values found to decrease with increase in material strength. Also, the flaw angle 45° was found to have the relatively high difference between these two stresses compared to other flaw angles. The numerical simulation could successfully predict the peak stress along after crack initiation stress and wing crack propagation along with subsidiary cracks similar to laboratory experiments.
Fig. 20 Comparison between the experimental and the numerical results of peak stress observed for three different strengths specimens with respect to flaw angle
References The difference between the values of the crack initiation stress and peak stress obtained from the experimental and numerical studies are quite higher for L30, L60, M30 and M60 compared to other flaw angles. This is due to the presence of both Type 1 and Type 2 wing crack which is obeserved for the same flaw angle specimen in the experiments. But in the numerical study for the same flaw angle either Type 1 or Type 2 or transformation between these
1. Griffith AA (1920) The phenomenon of rupture and flow in solids. Philos Trans R Soc Lond Ser A 221:163–198 2. Inglis CE (1913) Stresses in a plate due to the presence of cracks and sharp corners. Trans Inst Nav Archit 60:219–230 3. Orowan E (1949) Fracture and strength of solids. Rep Progr Phys 12:185–232 4. Hoek E (1964) Fracture of anisotropic rock. J S Afr Inst Min Metall 64(10):501–518
123
Indian Geotech J 5. Hoek E, Bieniawski ZT (1965) Brittle fracture propagation in rock under compression. J Fract Mech 1(3):139–155 6. Barenblatt GI (1962) The mathematical theory of equilibrium of cracks in brittle fracture. Adv Appl Mech 7:55–129 7. Dugdale DS (1960) Yielding of steel sheets containing slits. J Mech Phys Solids 8:100–104 8. Hillerborg A, Modeer M, Petersson PE (1976) Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 6:773–782 9. Bobet A, Einstein HH (1998) Numerical modeling of fracture coalescence in a model rock material. Int J Fract 92:221–252 10. Bobet A, Einstein HH (1998) Fracture coalescence in rock-type materials under uniaxial and biaxial compression. Int J Rock Mech Min Sci 35(7):863–888 11. Sagong M, Bobet A (2002) Coalescence of multiple flaws in a rock- model material in uniaxial compression. Int J Rock Mech Min Sci 39(2):229–241 12. Wong LNY, Einstein HH (2009) Crack coalescence in molded gypsum and Carrara marble: part 1—macroscopic observations and interpretation. Rock Mech Rock Eng 42(3):475–511 13. Wong LNY, Einstein HH (2009) Systematic evaluation of cracking behavior in specimens containing single flaws under uniaxial compression. Int J Rock Mech Min Sci 46(2):239–249 14. Park CH, Bobet A (2009) Crack coalescence in specimens with open and closed flaws: a comparison. Int J Rock Mech Min Sci 46(5):819–829 15. Park CH, Bobet A (2010) Crack initiation, propagation and coalescence from frictional flaws in uniaxial compression. Eng Fract Mech 77(14):2727–2748 16. Reyes O, Einstein HH (1991) Failure mechanism of fractured rock-a fracture coalescence model. In: Proceedings of the 7th congress of the ISRM, Aachen, Germany, vol 1, pp 333–340 17. Xu Y, Yuan H (2011) Applications of normal stress dominated cohesive zone models for mixed-mode crack simulation based on extended finite element methods. Eng Fract Mech 78:544–558 18. Gonc¸laves da Silva B, Einstein HH (2013) Modeling of crack initiation, propagation and coalescence in rocks. Int J Fract 182:167–186 19. Xie Y, Cao P, Liu J, Dong L (2016) Influence of crack surface friction on crack initiation and propagation: a numerical investigation based on extended finite element method. Comput Geotech 74:1–14 20. Lee H, Jeon S (2011) An experimental and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression. Int J Solids Struct 48(6):979–999 21. Lisjak A, Figi D, Grasselli G (2014) Fracture development around deep underground excavations: insights from FDEM modelling. J Rock Mech Geotech Eng 6:493–505 22. Mahabadi OK, Lisjak A, Grasselli G, Munjiza A (2012) Y-Geo: a new combined finite-discrete element numerical code for geomechanical applications. Int J Geomech 12(6):676–688
123
23. Vasarhelyi B, Bobet A (2000) Modeling of crack initiation, propagation and coalescence in uniaxial compression. Rock Mech Min Sci 33(2):119–139 24. Steif PS (1984) Crack extension under compressive loading. Eng Fract Mech 20(3):463–473 25. Horri H, Nemat-Naseer S (1985) Compression-induced microcrack growth in brittle soilds: axial splitting and shear failure. J Geophys Res 90(B4):3105–3125 26. Ashby MFA, Hallam SD (1986) The failure of brittle solids containing small cracks under compressive stress states. Acta Metall Sin 34(3):497–510 27. Baud P, Reuschle´ T, Charlez P (1996) An improved wing crack model for the deformation and failure of rock in compression. Int J Rock Mech Min Sci 33(5):539–542 28. Lajtai EZ (1974) Brittle fracture in compression. Int J Fract 10(4):525–536 29. Xu XP, Needleman A (1994) Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 42:1397–1434 30. Camanho PP, Da´vila CG (2002) Mixed-mode decohesion finite elements for the simulation of delamination in composite materials. NASA/TM-2002-211737, pp 1–42 31. Belytschko T, Black T (1999) Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 45:601–620 32. Fairhurst CE, Hudson JA (1999) Draft ISRM suggested method for the complete stress–strain curve for intact rock in uniaxial compression. Int J Rock Mech Min Sci 36(3):279–289 33. ISRM (1979) Suggested methods for determining the uniaxial compressive strength and deformability of rock materials. Int J Rock Mech Min Sci Geomech Abstr 16(2):135–140 34. ISRM (1978) Suggested methods for determining tensile strength of rock materials. Int J Rock Mech Sci Geomech Abstr 15:99–103 35. ISRM (2007) The complete ISRM suggested methods for rock characterization, testing and monitoring: 1974–2006. In: Ulusay R Hudson JA (eds) Suggested methods prepared by the commission on testing methods, ISRM, compilation arranged by the ISRM Turkish National Group, Kozan Ofset, Ankara 36. Fowell RJ (1995) Draft ISRM suggested method for determining mode I fracture toughness using cracked chevron notched barzilian disc (CCNBD) specimens. Int J Rock Mech Min Sci Geomech Abstr 32(1):57–64 37. Atkinson C, Smelser RE, Sanchez J (1982) Combined mode fracture via the cracked Brazilian disk test. Int J Fract 18(4):279–291 38. Al-Shayea NA (2005) Crack propagation trajectories for rocks under mixed mode I–II fracture. Eng Geol 81:84–97 39. Yamashita T, Umeda Y (1994) Earthquake rupture complexity due to dynamic nucleation and interaction of subsidiary faults. Pure Appl Geophys 143(1–3):89–116 40. Abaqus Analysis User’s Manual, Version 6.12 (2012) Dassault Syste`mes Simulia Corporation. Providence, Rhode Island, USA