- Email: [email protected]
Discrete element simulation of SiC ceramic containing a single pre-existing flaw under uniaxial compression
Author’s Accepted Manuscript Discrete element simulation of SiC ceramic containing a single pre-existing flaw under uniaxial compression Shengqiang Jiang, Xu Li, Li Zhang, Yuanqiang Tan, Ruitao Peng, Rui Chen www.elsevier.com/locate/ceri
PII: DOI: Reference:
S0272-8842(17)32556-7 http://dx.doi.org/10.1016/j.ceramint.2017.11.099 CERI16754
To appear in: Ceramics International Received date: 27 September 2017 Revised date: 14 November 2017 Accepted date: 14 November 2017 Cite this article as: Shengqiang Jiang, Xu Li, Li Zhang, Yuanqiang Tan, Ruitao Peng and Rui Chen, Discrete element simulation of SiC ceramic containing a single pre-existing flaw under uniaxial compression, Ceramics International, http://dx.doi.org/10.1016/j.ceramint.2017.11.099 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Discrete element simulation of SiC ceramic containing a single pre-existing flaw under uniaxial compression
Shengqiang Jiang a,*, Xu Lia, Li Zhanga, Yuanqiang Tanb, Ruitao Penga and Rui Chena a
School of Mechanical Engineering, Xiangtan University, Hunan, 411105, China
b
Institute of Manufacturing Engineering, Huaqiao University, Fujian, 361021, China
* Corresponding author. Tel.: +86 731 58292225. Fax: +86 731 58292209. E-mail address: [email protected]
Abstract: A model of a SiC ceramic containing a single pre-existing flaw was established based on the discrete element method. The effects of the flaw inclination angles, which ranged from 0° to 75°, on the mechanical properties of the specimen under uniaxial compression were studied. The evolution of the force-chain field, displacement field and stress field around the pre-existing flaw in the process from the load to failure was also analysed. The results showed that the flaw inclination angle affected the mechanical properties of the specimen as well as the initiation and propagation of the first crack. Based on the investigation of the force chain field, it was found that the distribution curve of the normal force carried by the parallel bond in the specimen with the corresponding angles under compression is similar to the “peanut” rose diagram, while the shear force distribution curve is similar to the "butterfly wings" rose diagram. In addition, in the analysis of the displacement field and the stress field, the displacement field around the flaw can be divided into four types in the process from specimen loading to its failure. Meanwhile, it was found that initiation of the first crack was affected by tensile stress. With the propagation of the first crack, the tensile stress concentration region at the flaw tip moved and dissipated correspondingly.
Keywords˖Crack propagation mechanism; Single pre-existing flaw; Uniaxial compression; Discrete element method
1. Introduction In traditional mechanics, materials are considered to be ideal solids without any defects. In fact, a variety of cracks are induced in the preparation, processing and service of materials [1]. As a result, the presence of cracks can lead to a certain degree of dispersion of the material strength and reduce the reliability of the material. Moreover, the defect sensitivity is more salient, and an unpredictable brittle fracture during service occurs easily in engineering ceramics and other hard and brittle materials. Therefore, for the study of failure mechanisms, it is important to explore the process of interval crack initiation, propagation and coalescence in hard and brittle materials such as the engineering ceramics
during their preparation, processing and service [2-4]. In previous studies, hard and brittle materials such as rock [5-9] and gypsum [10-12] were selected as the main object of investigation, and a lot of research has been focused on the cracking process of materials under pressure, with their results mainly based on experimental measurements and numerical simulations. In experimental studies, most researchers used model tests to study the mechanism of crack propagation through the use of acoustic emission (AE) [13] and computerized tomography (CT) [14] techniques to observe and record the cracking process. However, there is a considerable difference between the model sample and the actual material, which will inevitably lead to the inaccuracy of the study. Meanwhile, the interval flaws are mostly produced by machining and embedding, generating a large pre-existing flaw and inevitably inducing new damage. In this field, the discrete element method (DEM) as a numerical simulation method has attracted increasing research attention in recent years because of its unique advantages for the study of the crack of brittle materials. Proposed by Cundall [15] for the study of rock mechanics in the early 1970s, the DEM has been widely used in the study of crack initiation, propagation and coalescence of rock materials [16-18] and has been applied further to the study of the fracture, failure and processing of hard brittle materials such as concrete [19], laminated glass [20], and ceramics [21-23]. Our previous study [24] has shown that the presence of random microcracks in the ceramic material will affect its mechanical properties (such as uniaxial compressive strength, Young’s modulus and Poisson's ratio) and has a significant and decisive influence on the failure mode of the specimen. However, the failure of brittle materials, such as ceramics, is inseparable from the initiation and propagation of cracks. In this work, we will use the discrete element method to study the internal crack propagation mechanism of SiC ceramic materials under uniaxial compression. The effects of the pre-existing flaw, with inclination angles ranging from 0° to 75°, on the mechanical properties of the specimen were studied by the established DEM model of the SiC ceramic containing a single linear oblique flaw under uniaxial compression. At the same time, the evolution process of the force-chain field, displacement field and stress field around the pre-existing flaw in the process from loading to failure of the sample was also analysed.
2 . Methodology 2.1 Discrete element method Discrete element method is a numerical simulation method proposed by Cundall [15] in the early 1970s. It treats the material as a collection of discrete rigid particles which are connected by a bond. Potyondy and Cundall designed the bonded-particle model (BPM) [25] to simulate the mechanical behavior of bulk materials. In the BPM, it assumes that there exists a finite adhesive between two rigid particles, which can simultaneously transfer the force and moment between the particles. The assumed adhesive is represented by a parallel bond, and contacted particles are bonded together by means of parallel bonds to form an arbitrary shape assembly. The initiation and propagation of microcracks
in assembly can be seen as a breakage of the parallel bond between contacted particles. In this way, it can intuitively deal with the initiation and propagation process of crack, as shown in Fig. 1. In the DEM, the interaction between the contacted particles is a dynamic process, and the movement of particles at any time step meets Newton's second law. The movement has two ways of translation and rotation: Translation: Fi E g vi
Rotation: M i E g wi
m
I
'vi 't
'wi 't
(1)
(2)
where i (=1, 2) denotes the x and y coordinate directions, respectively; Fi is the out-of-balance force component of the particle; vi is the translational velocity; m is the mass of the particle; Mi is the out-of-balance moment due to the contacts; wi is the rotational velocity; I is the rotational inertia of the particle; βg is the global damping coefficient; Δt is the time step. Two arbitrary contacted particles are governed by the force-displacement rule, the contact force F is expressed as:
Fnn Fss
F
(3)
where Fn and Fs are the shear contact force and normal contact force respectively, the calculation formulas are as follows:
Fn 'F s
K nU n
(4)
K s 'U s
(5)
where Kn and Ks are normal and shear stiffness values at the contact, respectively; Un is the overlap of the two contacting particles; Fs is the increment of the shear contact force; ΔUs is the relative shear displacement (note that the shear contact force is initialised to zero when the two particles are in contact). The maximum tensile and shear stresses acting on the bond periphery are calculated (via beam theory) to be:
Vm a x
n M F R A I
(6)
s
Wm a x
Fi A
(7)
where A is the area of the cross section of a bond; I is the moment of inertia of the cross section of the bond; CF and M denote the total force and moment associated with the parallel bond respectively; CR is the radius of bond. The maximum tensile stress exceeds the tensile strength (σmax ≥ σbond) or the maximum shear stress exceeds the shear strength (τmax ≥ τbond), then the parallel bond breaks. And the relevant force, moment and stiffness are removed [27].
2.2 Modeling and calibration By PFC2D software developed by Itasca Company, two BPM models in which the particles were randomly distributed and closely arranged by adopting discrete element method was created in this paper. The dimensions of model are 4 × 8 mm and 2.4 × 0.5 mm, respectively, and the minimum radius of particles rlo was 1.5×10-2mm. Then, in order that the macroscopic mechanical properties of the model were consistent well with the actual SiC ceramic materials, as shown in Fig. 2, uniaxial compression test (the dimension of model is 4 × 8 mm), single edge notched beam test (the dimension is 2.4 × 0.5 mm, the span is 2 mm, the height of incision is 0.25 mm) and three-point bending test (the dimension is 2.4 × 0.5 mm, the span is 2 mm) were simulated respectively. By the trial and error method, the microscopic parameters of the discrete element model were constantly debugged and calibrated. The calibration results are shown in Table 1, it can be seen that the DEM simulation results are in good agreement with the experimental results. The discrete element model parameters of SiC ceramic are shown in Table 2. Table 1 Main mechanical properties of SiC in experimental tests and in DEM simulations Mechanical properties
Experimental results[21]
DEM simulation
Elastic Modulus (GPa)
420
422
0.14
0.14
Compressive strength (MPa)
2000
1935
Bending strength (MPa)
500-800
671
Fracture toughness (MPa/m1/2)
3.5
3.7
Poisson's ratio
Table 2 Parameters in the DEM model of SiC ceramic Variable name(unit)
Symbol
Value
Minimum radius of particles (m)
rlo
1.5e-5
Maximum-to-minimum radius ratio of the particle
n
1.2
Density of particle (kg/m3)
ρ
2600
Effective modulus of particle(Pa)
Eball
234e9
Normal-to-shear stiffness ratio of particle
kball
1.22
Effective modulus of parallel bond (Pa)
Ebond
234e9
Normal-to-shear stiffness ratio of parallel bond
kbond
1.22
Tensile strength of parallel bond (Pa)
σbond
9.3e8
Shear strength of parallel bond (Pa)
τbond
78e8
Friction coefficient of particles
μ
0.5
Deleted particles number of the pre-existing flaw
n
26
2.3 DEM model containing a single pre-existing flaw In order to produce a prefabricated central oblique flaw in the DEM model of SiC ceramic and ensure the width of resulting pre-existing flaw is as small as possible to meet the crack width distribution of the actual material at the same time, the method of deleting the particles was chosen. The diameter of the deleted particles represented the width of the flaw. And the half-length of the flaw can be obtained from:
a nu
d 2
n u rlo
(8)
where n is the deleted particles number along with flaw direction (as shown in Table 2), d is the minimum diameter of the particle, rlo is the minimum radius of the particle. In order to center the central oblique flaw with the center of the model, the length of the central oblique flaw can be determined by the deleted particles number along with flaw direction. The inclination angle θ of the center oblique flaw was represented by the included angle between the flaw direction and the positive direction of the x-axis. Moreover, θ can be arbitrarily specified as a controllable variable in the subsequent simulation. The model diagram is shown in Fig. 3.
2.4 Uniaxial compression test A uniaxial compression test of the discrete element model of SiC ceramics containing a single pre-existing flaw was established, as shown in Fig. 4. The pressure in the Y direction was applied to the specimen by setting a wall at the top and bottom of the discrete element model. The wall was represented by a horizontal line, whose stiffness was much larger than those of particles to prevent its deformation during the loading process. During the loading, both the top and bottom walls were moved to the center of the specimen with a constant v = 0.2 m/s by servo control. When the measured specimen stress is lower than 80% of the peak stress (defined as the compressive strength of the material), the loading was stopped and the specimen was broken finally. The inclination angle θ of the center oblique flaw in the specimen ranged from 0° to 75°. With 5° increment of inclination angle, those angles can be divided into 16 groups. So the fracture picture of the DEM model at different inclination angles was shown in Fig. 5. In those figures, the white straight line is the pre-existing center oblique flaw, while the red part indicates the microcracks generated in the loading process.
3. Results and Discussion
3.1 Analysis of Crack propagation During the loading process, the specimen went through four stages: initial compression, initiation of the first crack, generation of the secondary crack and specimen failure; for convenience, “stage 1”, “stage 2”, “stage 3” and “stage 4”, respectively, are used to represent these four stages. In Fig. 5 (stage 4), it was observed that the macroscopic crack angle was roughly between 60° and 70° when the specimen was broken at different flaw inclination angles. The failure modes of the specimen were basically the same. The break surface of the specimen that generated the macroscopic cracks produced many microcracks and resulted in many exfoliated blocks consisting of particles and exfoliated particles, as shown in Fig. 6. It is easy to infer that the largest stress is exerted on the specimen in the 60°~70° angle range. At the same time, there was a shear slip in this interval that finally resulted in a destabilizing rupture. Fig. 7 shows the initiation position of the first crack and secondary crack around the pre-existing flaw under compression. To eliminate the random error caused by the particle arrangement, we carried out six group numerical simulations on the specimen with the flaw inclination angles in the range of 0°~75° (interval of 5°). The enlarged images of the first crack and the secondary crack were chosen with the flaw angles at 0°, 25°, 45° and 65°, respectively, as shown in Figs. 8 and 9. It can be observed that the initiation position of the first crack was slightly different at the same flaw angle in Fig. 8. Additionally, the two growth tracks of the first crack were not symmetrical. With the increase of the flaw inclination angle θ, the distances d1 and d2 (the distance between the point of first crack initiation and flaw tip, respectively) were gradually reduced. At the same time, the two initiation angles of the first crack α1 and α2 increased (from approximately 90° to approximately 130°). These findings are in good agreement with the results of previous experiments [28] and simulations [16, 18]. Meanwhile, the initiation and propagation of the first crack were stable. A greater angle of the first crack makes it more difficult for the first crack to expand, and the propagation length of the crack was also gradually reduced. It has been widely confirmed in experimental studies of various brittle materials [30-33] that the initiation of the secondary crack at the flaw tips is normally subjected to shear force. Meanwhile, the secondary crack includes a coplanar secondary crack and an anti-wing crack, as shown in Fig. 7. The corresponding directions of the cracks were coplanar or nearly coplanar to the flaw and with an inclination similar to the wing cracks but in the opposite direction. Fig. 9 shows the enlarged images of the secondary crack when the flaw inclination angles are 0°, 25°, 45° and 65° (six distribution simulations were made for each flaw inclination angle.) It can be seen from the figure that the initiation position and direction of the secondary crack are in good agreement with the previous literature. The initiation and propagation of the secondary crack were unstable, and its crack surfaces were rough with some particles exfoliated by shearing.
3.2 Analysis of mechanical properties
The pre-existing flaw affected the mechanical properties of the specimen to a certain extent. Fig. 10 shows the DEM simulation results of the mechanical properties of SiC ceramics with a flaw for which the inclination angle ranged from 0° to 75° (the results were the average of the results of six different simulations, and the error bars represents the relevant standard deviation). The horizontal and vertical axes represent the flaw inclination angle and the ratio of the crack initiation stress to the compressive strength of the specimen, respectively. By comparing to the experimental results of for gypsum [28] and the DEM simulation results for a rock-like material [17], it can be concluded that the corresponding trend was roughly the same. At the same time, the ratios in the rock-like material and the SiC ceramics in the DEM simulation results were smaller than that obtained experimentally for the Gypsum material. One plausible explanation is that this was due to experimental limitations. In the physical tests, the detection of the first crack was based on the first distinct cracking sound or the minute crack opening event observed on the specimen surface. As suggested in other experimental studies employing the AE and CT techniques, cracking events prior to those observed on the surface could have occurred inside the specimens. A proper identification of the latter could lead to the lower ratio [17]. Fig. 11 shows the trends in Young's modulus and Poisson's ratio of the specimen when the flaw inclination angle ranged from 0° to 75°. It can be seen that the presence of the flaw in the specimen reduces the effective Young's modulus (420 GPa with no flaw) of the specimen and increases the Poisson's ratio (0.14 with no flaw) obtained in the uniaxial compression test. With the increase of the flaw inclination angle, the effective Young's modulus changed from being relatively stable to increasing, but the obtained value was always smaller than the original value. Meanwhile, Poisson's ratio was first relatively stable and then decreased, but it always remained higher than its original value.
3.3 Analysis of force-chain Fig. 12 shows the distribution of the force-chain field around the flaw at four different stages for the flaw inclination angles of 0°, 25°, 45° and 65° (stages 1, 2, 3 and 4 are the processes of initial compression, initiation of the first crack, generation of a secondary crack and specimen failure, respectively). In the figure, the particles are hidden, and the green and the black lines represent the tensile and the compressive force-chains of the parallel bond, respectively. Additionally, the thickness of the force-chain lines represents the magnitude of the contact resultant force on the parallel bond, while its direction indicated the direction of the contact resultant force. Finally, the red represented the microcracks generated during the loading process. Because in the DEM simulation the generation of microcracks is characterized by the breaking of the parallel bonds between the particles, the magnitude and direction of the force on the parallel bond around the flaw in the loading process (the distribution of the force-chain field) can, to a certain extent, reflect the force distribution and evolution around the flaw. Fig. 13 shows the evolution process of the strong force-chain field around the flaw for the inclination angles of 0°,
25°, 45° and 65° at different loading stages. Combining Figs. 12 and 13, we can observe that at the initial compression stage (stage 1), the concentration regions of the tensile force-chain (the green region A) were generated in the top and bottom of the pre-existing flaw at different inclination angles. Meanwhile, the concentration regions of the compressive force-chain (the black region B) were generated in the left and right edges of the pre-existing flaw. This finding indicated that a stress concentration was generated in regions A and B. However, the position and size of the concentration regions of the tensile force-chain (the green region A) changed for different flaw inclination angles. For example, for the flaw inclination angle of 0°, the A regions were just above or below the flaw, and the A region area was the largest. With the increase of the flaw inclination angle, the regions were gradually reduced and transferred to the flaw tip. By comparing the difference of the first crack initiation position in section 3.1 above, we can draw the conclusion that the first crack initiation position occurred in the concentration regions of the tensile force-chain (the green region A). With the decrease of the flaw inclination angle, both the concentration regions of the tensile force-chain and the distribution interval of the first crack initiation position become larger. Thus, the decrease of the flaw inclination angle in the experiments reported in the literature [29] increased the distribution error of the distance d1 and d2 between the first crack initiation position and the flaw tip instead. With the constant loading, the first crack in stage 2 was generated in the concentration regions of the tensile force-chain (the green region A) and then was consistently extended. Concentrated tensile stress was released during the first crack initiation and extension. Meanwhile, the region A transferred to the first crack tip and gradually shrank. When the secondary crack occurred in stage 3, the first crack stopped extending, and the tensile force-chain concentration region (the green region A) disappeared. However, this process did not have a significant effect on the compressive force-chain concentration region (the black region B) at both ends of the flaw. In addition, around the flaw, a new region called the “force-chain dissipation region” (the blue region C) was generated by the initiation of the first crack and changed accordingly with the extension of the first crack. It can be seen from Fig. 13 that the increasing flaw inclination angle made it difficult for the first crack to extend. As a result, the length of the first crack length decreased and the area of region C decreased as well. To a certain extent, the “force-chain dissipation region” (the blue region C) ensured the steady extension of the first crack as well as the absence of a new crack branch in the process of first crack extension. A quantitative analysis of the force-chain in the specimen was conducted in the following figures. Fig. 14 shows the distribution between the normal force magnitude of the parallel bond and the force-chain angle at different loading stages when the flaw inclination angles were 0°, 25°, 45° and 65°. The positive value of the normal force denoted compression while the negative value denoted tension. The force-chain angle was the included angle between the normal unit vector direction of the parallel bond and the positive direction of the x-axis. In Figs. 14(a-d), in the initial stage, the normal force of the parallel bond was initialized to zero. As the loading was carried out, the normal force
began to increase, and the top and bottom of the corresponding curve also began to extend in the 90° and 270° directions, indicating the maximums of the normal compression force on the parallel bond around the angle intervals at that time. However, the left and right ends of the curve were retracted in the opposite direction towards 180° and 0°, indicating that in the angle intervals, the normal tension force on the parallel bond was maximized. Finally, the final curve gradually became a “peanut” rose diagram. When the specimen failed in stage 4, the value of the normal force on the parallel bond showed a drop at the different flaw inclination angles, and the curve moved back by a certain amount. Meanwhile, it also found that the curve was sunken slightly in the 90° and 270° direction, and the normal force on the parallel bond decreased slightly in the specimen. It is possible that the propagation and coalescence of the crack in the specimen weakens the load-bearing capacity of the specimen at the angle intervals. Fig. 15 shows the distribution between the shear force magnitude of the parallel bond and the force-chain angles at the different loading stages in the specimen for the flaw inclination angles of 0°, 25°, 45° and 65°. As seen from Figs. 15(a-d), the shear force on the parallel bond is initialized to zero. As the load increased, the shear force increased accordingly. The curves in the angle intervals in the 50°-55°, 140°-145°, 230°-235° and 320°-325° intervals continue to extend, and the shear force reached maximum values in the four angle intervals. As a result, the final curve gradually became the "butterfly wings" rose diagram. When the specimen failed in stage 4, the shear force magnitudes of the parallel bond dropped at the different flaw inclination angles, and the curve moved back by a certain amount back as well. The comparison between the above study and the angle interval range of the generated macrocracks in the previous section 3.1 confirmed the internal origin of the macroscopic cracks generated in the 60° ~ 70° angle interval when the specimen failed.
3.4 Analysis of displacement During the process from initial compression to rupture failure of the discrete element specimen, the displacement of the particles can be indicated by hiding the particles and only showing the displacement tendency using arrows. To some extent, this facilitated the observation of the process of microcrack initiation, propagation and coalescence at different loading stages. The displacement tendency arrows show a corresponding evolution at different loading stages. As a result, we divided the displacement field into four different types, denoted type I, type II, type III, and type IV for convenience. As shown in Fig. 16, type I represents a purely tensile displacement field, type II represents both the tensile and shear displacement field, type III represents a displacement field that was neither tensile nor shear, and type IV represents both the tensile and shear displacement fields but with the tensile field much larger than the shear field. Figs. 17 (a-d) show the evolution process of the displacement tendency arrows for the flaw inclination angles of 0°, 25°, 45° and 65°. It can be seen that the arrows around the first crack were open to both sides of the first crack regardless of the change of the flaw inclination angles at stage 2. At that time, the displacement field of the first crack
was type Ϩ, indicating that the first crack was generated by tension. As the loading progressed, the secondary cracks in stage 3 were generated accordingly. The displacement tendency arrows around the crack initiation position gradually moved from their original upward or downward positions towards the direction far away from both sides of flaw tips after the arrows approached the flaw tips. The displacement field of the secondary cracks was almost identical to type II, which meant that the secondary crack was generated by both tension and shear. At the same time, we also observed that the displacement tendency arrows around the first crack began to evolve from type Ϩ to type Ϫ because of secondary crack generation and some of first cracks gradually stopped extending. During the constant loading process, some of the secondary cracks gradually developed into macroscopic cracks. When the specimen was completely broken in the stage 4, the displacement field showed more changes. The changes were as follows: the displacement tendency arrows around the first crack almost turned into type III because the first crack no longer expanded while the arrows around the secondary crack shifted from the original type ϩ to type Ϫor type ϫ. The displacement field where the secondary crack evolved into the macroscopic crack changed from the original type II to type ϫ, and the displacement field where the secondary crack stopped the initiation transformed into type III.
3.5 Analysis of stress The study of the evolution of the stress field around the flaw at different loading stages is beneficial for further understanding of the entire process of crack initiation, propagation and convergence. By setting up a series of measurement circles around the flaw, the measurement of the stress distribution around the flaw was carried out in this work and then a stress distribution cloud picture was drawn. The measured stress magnitude was the average of the stress carried by the particles contained in the circle. The number of particles was generally 4-6. The distribution of the measurement circle is shown in Fig. 18. Here, the stress distributions in stages 1, 2 and 3 were selected and measured for the crack inclination angles of 0°, 25°, 45° and 65°. Then, the stress distribution cloud picture of the three stress components, namely, σxx, σxy and σyy, were generated accordingly. As shown in Figs. 19-21, the coordinate system in the stress cloud picture adopted a local coordinate system of the flaw. The centre of the flaw was the origin of the coordinate system, and the x- and y-coordinate values divided by the half-length of the preset crack were the dimensionless numbers. A negative stress value denoted compression and a positive value corresponded to tension, in the units of MPa. Figs. 19(a-d) showed the σxx stress component distribution around the flaw at different stage for the flaw inclination angles of 0°, 25°, 45° and 65°. For the flaw inclination angle of 0°, we observed that the compressive stress concentration regions (blue regions) appeared at the flaw tips. Meanwhile, a large range of tensile stress concentration regions (red regions) appeared immediately above and below the flaw. With the increase of the flaw inclination angle, the compressive stress concentration regions at the flaw tips moved correspondingly to the centre of the flaw.
Meanwhile, the tensile stress concentration regions above and below the flaw also moved towards the flaw tips. In the process of changing the inclination angle of the flaw from 0° to 65°, we found that both the compressive stress concentration regions and the tensile stress concentration regions around the flaw exchanged their positions. With the loading, the first crack was generated and extended in the stage 2. It can be observed in the figure that the first crack initiated in the tensile stress concentration regions (red regions) and released the tensile stress in this region. With the extension of the first crack, the original tensile stress concentration regions transferred to the first crack tips, while the compressive stress concentration regions (blue regions) were not affected. Then, the secondary crack initiated and extended near the flaw tips in stage 3, accelerating the dissipation of the tensile stress concentration regions. However, there was no effect on the tensile stress concentration regions. Fig. 20 shows the σxy stress component distribution around the flaw with different inclination angles at different stages. The transfer phenomenon of the stress concentration regions was also observed in stage 1 when the specimen was loaded at the same time. In the process of changing the inclination angle of the flaw from 0° to 65°, the positive stress concentration regions were initially located in the left bottom and right top of the flaw rotated clockwise 180° around the flaw tip and then transferred to the left top and right bottom of the flaw. The negative stress concentration region was initially located in the left top and right bottom of the flaw and rotated clockwise around the flaw and transferred to the right top and left bottom of the flaw. It is important to note that the initiation and expansion of the first and secondary cracks in stages 2 and 3 did not affect the evolution of the σxy stress component. Fig. 21 shows the distribution of the stress component σyy around the corresponding flaw. The stress component σyy in the stress cloud picture was negative due to the compression loading of the specimen. It was obvious that the compressive stress concentration regions were generated at the flaw tips (blue regions). Meanwhile, it appears shielding regions of compressive stress (brown regions) were formed at the top and bottom of the flaw. The shielding effect was stronger at smaller distances from the flaw surface. As a result, the compressive stress became weaker.
3.6 Discussion and outlook As shown in Figs. 10 and 11, it is important to note that when the flaw inclination angle was in the range of 30-40º, all curves of the ratio of crack initiation stress to compressive strength, Young's modulus and Poisson's ratio changed. The possible reasons for the change are as follows: (1) for the length of first crack, when the flaw inclination angle was in the range of 0-30º, the first crack was easier to extend. This characteristic resulted in the greater length of the first crack when the specimen failed under loading (as shown in Figs. 8 and 9). However, for the flaw inclination angles greater than 30º, it became more difficult for the first crack to extend. Therefore, the length of the first crack for the flaw inclination angles greater than 30º was relatively shorter. The difference in the length of the first crack resulted in the difference in the energy consumption at the flaw tips. Accordingly, the mechanical properties of the specimen were
different. (2) For the size and dissipation of the concentration regions of the tensile force-chain and tensile stress, when the flaw inclination angle was in the range of 0-30º, the concentration regions of the tensile force-chain (green regions in Figs. 12 and 13) and the tensile stress (red regions in Fig. 19) near the flaw tips were larger. The tensile stress concentration phenomenon at the flaw tips was more distinct than for the cases where the flaw inclination angle was greater than 30º. As a result, when the first crack was generated and extended, the “force-chain dissipation region” was larger, and the tensile stress dissipation effect was more obvious. In this work, a model of a SiC ceramic containing a single pre-existing flaw was established based on the discrete element method. However, the distribution of the microcrack is random and their direction cannot be controlled in the process of preparation, processing and service of the actual material. For engineering ceramics, the microcracks caused by machining damage and intrinsic defects can be detected and measured by non-destructive testing techniques (penetrant testing, ultrasonic testing, laser ultrasonication and infrared testing). The different flaw inclination angles will result in the differences of the mechanical properties of the specimens. To some extent, we can predict and evaluate the mechanical properties of the specimen using the DEM modelling based on the reasonable statistics for the distribution of the microcrack with different inclination angles in the specimen. Thus, this work can increase the reliability and service life of the specimen under different service conditions. At the same time, ceramic materials show a strong sensitivity to crack defects and a high fragility. An unpredictable brittle fracture in service is mostly caused by the initiation, propagation and coalescence of interval cracks. Ceramic toughening techniques can effectively improve the fracture toughness of ceramics as well as their reliability and service life. To enhance the toughness of ceramic materials, it is crucial to improve their ability to resist crack propagation and to reduce the stress concentration effect of crack tips. Studies of the relationship between the flaw inclination angle and the location and angle orientation of the toughening material in specimen can further reveal the inherent mechanism of ceramic toughening and explore the effect of the location and arrangement of toughening materials around the flaw on the toughening effect of the specimen under the conditions of different flaw inclination angles. Additionally, such studies can provide guidance for the efficient toughening of ceramic materials.
4
Conclusion (1) The flaw reduces the effective Young's modulus of the specimen and increases the Poisson's ratio under the
uniaxial compression test. With the increase of the flaw inclination angle, the effective Young's modulus changes from a being relatively stable to increasing, but its value is always smaller than the original value. Meanwhile, Poisson's ratio is first relatively stable and then decreases, and this value is always greater than its original value. (2) In the process of changing the inclination angle of the flaw from 0° to 65°, the concentration regions of the tensile force-chain above and below the flaw are gradually reduced and transferred to the flaw tip. The first crack
generated in the concentration regions of the tensile force-chain makes the region transfer to the first crack tip and decrease gradually. (3) During the loading process, the distribution curve between the normal force magnitude of the parallel bond and the force-chain angle gradually transforms into the "peanut" rose diagram. The shear force curve with angles gradually becomes similar to the "butterfly wings" rose diagram. When the specimen fails, the force magnitude of the parallel bond shows a drop. The study of the normal and shear forces of the parallel bond with the angle distribution in the specimen confirms the intrinsic origin of the macroscopic crack generation in the 60°-70° angle interval. (4) The evolution of the displacement field around the flaw in the process from the first crack initiation to the specimen failure is analysed. The initiation, propagation and coalescence of the crack are summarized by four types of displacement tendency arrows (type I, type II, type III and type IV). It is confirmed that the initiation of the first crack is affected by the tensile force while the secondary crack is affected by both the tensile and shear forces. (5) The first crack initiated in the tensile stress concentration regions. With the extension of the first crack, the tensile stress in the region is released. The tensile stress concentration regions transfer to the first crack tip and decrease. However, the initiation and propagation of the first and secondary cracks have little effect on the compressive stress concentration regions at the flaw tips. At the same time, examination of stress cloud picture analysis of the σxx and σxy stress components around the flaw under initial compression shows that the positive and negative stress concentration regions are transferred clockwise along the flaw in the process of changing the flaw inclination angles from 0° to 65°. In the stress cloud picture analysis of the σyy stress component, the compressive stress concentration is generated at the flaw tips. It appeared to form shielding regions of compressive stress at the top and bottom of the flaw. Stronger shielding effects were observed closer to the flaw surface. As a result, the compressive stress became weaker.
Acknowledgement This work was supported by the National Natural Science Foundation of China (Grant No. 51605409, 11772135, 51475404 and 11602212)
References [1] B. Lawn, Fracture of brittle solids, Cambridge university press, Cambridge, 1993. [2] H. Abe, M. Naito, T. Hotta, Flaw size distribution in high-quality alumina, J. Am. Ceram. Soc. 86 (2010) 1019-1021. [3] W.H. Tuan, J.C. Kuo, Effect of abrasive grinding on the strength and reliability of alumina, J. Eur. Ceram. Soc. 18 (1998) 799-806. [4] K. Ikeda, H. Igaki, Effect of surface flaw size on fracture strength of alumina ceramics, J. Am. Ceram. Soc. 70
(1987) C-29-C-30. [5] Z.T. Bieniawski, Mechanism of brittle fracture of rock: parts 1, Int. J. Rock. Mech. Min. 4 (1967) 395-404. [6] L.N. Germanovich, R.L. Salganik, A.V. Dyskin, Mechanisms of brittle fracture of rock with pre-existing cracks in compression, Pure. Appl. Geophys. 143 (1994) 117-149. [7] A. Bobet, H.H. Einstein, Fracture coalescence in rock-type materials under uniaxial and biaxial compression, Int. J. Rock. Mech. Min. 35 (1998) 863-888. [8] P. Cao, T. Liu, C. Pu, Crack propagation and coalescence of brittle rock-like specimens with pre-existing cracks in compression, Eng. Geol. 187 (2015) 113-121. [9] C.H. Park, A. Bobet, Crack coalescence in specimens with open and closed flaws: a comparison, Int. J. Rock. Mech. Min. 46 (2009) 819-829. [10] B. Shen, The mechanism of fracture coalescence in compression-experimental study and numerical simulation, Eng. Fract. Mech. 51 (1995) 73-85. [11] L.N.Y. Wong, H.H. Einstein, Crack coalescence in molded gypsum and carrara marble: part 1-macroscopic observations and interpretation, Rock. Mech. Rock. Eng. 42 (2008) 475-511. [12] L.N.Y. Wong, H.H. Einstein, Crack coalescence in molded gypsum and carrara marble: part 2-microscopic observations and interpretation, Rock. Mech. Rock. Eng. 42 (2009) 513-545. [13] T.J. Ahrens, Rock physics & phase relations: A handbook of physical constants, American Geophysical Union, Washington DC, 1995. [14] X.T. Feng, S. Chen, H. Zhou, Real-time computerized tomography (CT) experiments on sandstone damage evolution during triaxial compression with chemical corrosion, Int. J. Rock. Mech. Min. 41 (2004) 181-192. [15] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Géotechnique. 29 (1979) 47-65. [16] X.P. Zhang, L.N.Y. Wong, Crack initiation, propagation and coalescence in rock-like material containing two flaws: a numerical study based on bonded-particle model approach, Rock. Mech. Rock. Eng. 46 (2013) 1001-1021. [17] X.P. Zhang, Cracking processes in rock-like material containing a single flaw under uniaxial compression: a numerical study based on parallel bonded-particle model approach, Rock. Mech. Rock. Eng, 45 (2012) 711-737. [18] L.A.M. Camones, E. do Amaral Vargas, R.P.D. Figueiredo, Application of the discrete element method for modeling of rock crack propagation and coalescence in the step-path failure mechanism, Eng. Geol. 153 (2013) 80-94. [19] S. Hentz, F.V. Donzé, L. Daudeville, Discrete element modeling of concrete submitted to dynamic loading at high strain rates, Comput. Struct. 82 (2004) 2509-2524. [20] M.Y. Zang, Z. Lei, S.F. Wang, Investigation of impact fracture behavior of automobile laminated glass by 3D discrete element method, Comput. Mech. 41 (2007) 73-83.
[21] Y.Q. Tan, D.M. Yang, Y. Sheng, Discrete element method (DEM) modeling of fracture and damage in machining process of polycrystalline SiC, J. Eur. Ceram. Soc. 29 (2009) 1029-1037. [22] Y.Q. Tan, D.M. Yang, Y. Sheng, Study of polycrystalline Al2O3 machining cracks using discrete element method, Int. J. Mach. Tool. Manu. 48 (2008) 975-982. [23] Y.Q. Tan, S.Q. Jiang, D.M. Yang, Scratching of Al2O3, under pre-stressing, J. Mater. Process. Tech. 211 (2011) 1217-1223. [24] S.Q. Jiang, X. Li, Y.Q. Tan, Discrete element simulation of SiC ceramic with pre-existing random flaws under uniaxial compression, Ceram. Int. 43 (2017) 13717-13728. [25] D.O. Potyondy, P.A. Cundall, A bonded-particle model for rock, Int. J. Rock. Mech. Min, 41 (2004) 1329-1364. [26] N. Cho, C.D. Martin, D.C. Sego, A clumped particle model for rock, Int. J. Rock. Mech. Min. 44 (2007) 997-1010. [27] P.A. Cundall, PFC2D user's manual (Version 3.1), Minnesota: Itasca Consulting Group Inc, Minnesota, 2004. [28] L.N.Y. Wong, H. Einstein, Fracturing behavior of prismatic specimens containing single flaws, Golden Rocks 2006, The 41st US Symposium on Rock Mechanics (USRMS), American Rock Mechanics Association, 2006. [29] A. Bobet, The initiation of secondary cracks in compression, Eng. Fract. Mech. 66 (2000) 187-219. [30] E.Z. Lajtai, Brittle fracture in compression, Int. J. Fracture. 10 (1974) 525-536. [31] A.R. Ingraffea, F.E. Heuze, Finite element models for rock fracture mechanics, Int. J. Numer. Anal. Met. 4 (1980) 25-43. [32] H. Jiefan, C. Ganglin, Z. Yonghong, An experimental study of the strain field development prior to failure of a marble plate under compression, Tectonophysics. 175 (1990) 269-284. [33] Y.P. Li, L.Z. Chen, Y.H. Wang, Experimental research on pre-cracked marble under compression, Int. J. Solids. Struct. 42 (2005) 2505-2516.
List of figure captions: Fig. 1. BPM model diagram in Ref.[26]. Fig. 2. Calibration tests of mechanical properties, (a) uniaxial compression test, (b) single edge notched beam test, and (c) three-point bending test. Fig. 3. (a) Model with a single flaw, (b) DEM model containing a single pre-existing flaw. Fig. 4. DEM model in uniaxial compression test. Fig. 5. Rupture picture of specimens at different flaw inclination angles (0°~75°). Fig. 6. (a) Macroscopic crack in specimen with 45° pre-existing flaw, (b) enlarged image of yellow circle “A” and (c) enlarged image of yellow circle “B”. Fig. 7. Initiation position diagram of the first crack and secondary crack (θ is the inclination angle of pre-existing flaw. α1 and α2 are the crack initiation angles. d1 and d2 are the distances between the point of first crack initiation and flaw tip)[28,29]. Fig. 8. Enlarged image of the first crack under different flaw inclination angles, (a) 0°, (b)25°, (c)45° and (d)65°. Fig. 9. Enlarged image of the secondary crack under different flaw inclination angle, (a) 0°, (b) 25°, (c) 45° and (d) 65°. Fig. 10. Comparison of the ratio of the crack initiation stress and compressive strength under different flaw inclination angles. Fig. 11. Young’s modulus and Poisson’s ratio under the different flaw inclination angles. Fig. 12. Distribution of the force-chain field around the flaw at four stages under the different flaw inclination angles (a) 0°, (b) 25°, (c) 45° and (d) 65° (green lines indicate the tensile force-chain, black lines indicate the compressive force-chain, and red lines indicate the ruptured parallel bonds). Fig. 13. Distribution evolution of the strong force-chain field around the flaw at different loading stages when the flaw inclination angles are (a) 0°, (b) 25°, (c) 45° and (d) 65° (green regions A indicate the concentration regions of the strong tensile force-chain, black regions B indicate the concentration regions of the strong compressive force-chain, and blue regions C indicate the “force-chain dissipation region”). Fig. 14. Contrast analysis of the normal force of parallel bond vs the angle under different flaw inclination angles. Fig. 15. Contrast analysis of the shear force of parallel bond vs the angle under different flaw inclination angles. Fig. 16. Four types of displacement tendency arrows at different loading stages. Fig. 17. Evolution process of displacement tendency arrows at different loading stages with the flaw inclination angles (a) 0°, (b) 25°, (c) 45°, and (d) 65°. Fig. 18. Measurement circle distribution diagram around the pre-existing flaw. Fig. 19. Stress component σxx around the flaw at different loading stages with different inclination angles: (a) 0°, (b) 25°, (c) 45°, and (d) 65°.
Fig. 20. Stress component σxy around the flaw at different loading stages with different inclination angles: (a) 0°, (b) 25°, (c) 45°, and (d) 65°. Fig. 21. Stress component σyy around the flaw at different loading stages with different inclination angles: (a) 0°, (b) 25°, (c) 45°, and (d) 65°.
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
PII: DOI: Reference:
S0272-8842(17)32556-7 http://dx.doi.org/10.1016/j.ceramint.2017.11.099 CERI16754
To appear in: Ceramics International Received date: 27 September 2017 Revised date: 14 November 2017 Accepted date: 14 November 2017 Cite this article as: Shengqiang Jiang, Xu Li, Li Zhang, Yuanqiang Tan, Ruitao Peng and Rui Chen, Discrete element simulation of SiC ceramic containing a single pre-existing flaw under uniaxial compression, Ceramics International, http://dx.doi.org/10.1016/j.ceramint.2017.11.099 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Discrete element simulation of SiC ceramic containing a single pre-existing flaw under uniaxial compression
Shengqiang Jiang a,*, Xu Lia, Li Zhanga, Yuanqiang Tanb, Ruitao Penga and Rui Chena a
School of Mechanical Engineering, Xiangtan University, Hunan, 411105, China
b
Institute of Manufacturing Engineering, Huaqiao University, Fujian, 361021, China
* Corresponding author. Tel.: +86 731 58292225. Fax: +86 731 58292209. E-mail address: [email protected]
Abstract: A model of a SiC ceramic containing a single pre-existing flaw was established based on the discrete element method. The effects of the flaw inclination angles, which ranged from 0° to 75°, on the mechanical properties of the specimen under uniaxial compression were studied. The evolution of the force-chain field, displacement field and stress field around the pre-existing flaw in the process from the load to failure was also analysed. The results showed that the flaw inclination angle affected the mechanical properties of the specimen as well as the initiation and propagation of the first crack. Based on the investigation of the force chain field, it was found that the distribution curve of the normal force carried by the parallel bond in the specimen with the corresponding angles under compression is similar to the “peanut” rose diagram, while the shear force distribution curve is similar to the "butterfly wings" rose diagram. In addition, in the analysis of the displacement field and the stress field, the displacement field around the flaw can be divided into four types in the process from specimen loading to its failure. Meanwhile, it was found that initiation of the first crack was affected by tensile stress. With the propagation of the first crack, the tensile stress concentration region at the flaw tip moved and dissipated correspondingly.
Keywords˖Crack propagation mechanism; Single pre-existing flaw; Uniaxial compression; Discrete element method
1. Introduction In traditional mechanics, materials are considered to be ideal solids without any defects. In fact, a variety of cracks are induced in the preparation, processing and service of materials [1]. As a result, the presence of cracks can lead to a certain degree of dispersion of the material strength and reduce the reliability of the material. Moreover, the defect sensitivity is more salient, and an unpredictable brittle fracture during service occurs easily in engineering ceramics and other hard and brittle materials. Therefore, for the study of failure mechanisms, it is important to explore the process of interval crack initiation, propagation and coalescence in hard and brittle materials such as the engineering ceramics
during their preparation, processing and service [2-4]. In previous studies, hard and brittle materials such as rock [5-9] and gypsum [10-12] were selected as the main object of investigation, and a lot of research has been focused on the cracking process of materials under pressure, with their results mainly based on experimental measurements and numerical simulations. In experimental studies, most researchers used model tests to study the mechanism of crack propagation through the use of acoustic emission (AE) [13] and computerized tomography (CT) [14] techniques to observe and record the cracking process. However, there is a considerable difference between the model sample and the actual material, which will inevitably lead to the inaccuracy of the study. Meanwhile, the interval flaws are mostly produced by machining and embedding, generating a large pre-existing flaw and inevitably inducing new damage. In this field, the discrete element method (DEM) as a numerical simulation method has attracted increasing research attention in recent years because of its unique advantages for the study of the crack of brittle materials. Proposed by Cundall [15] for the study of rock mechanics in the early 1970s, the DEM has been widely used in the study of crack initiation, propagation and coalescence of rock materials [16-18] and has been applied further to the study of the fracture, failure and processing of hard brittle materials such as concrete [19], laminated glass [20], and ceramics [21-23]. Our previous study [24] has shown that the presence of random microcracks in the ceramic material will affect its mechanical properties (such as uniaxial compressive strength, Young’s modulus and Poisson's ratio) and has a significant and decisive influence on the failure mode of the specimen. However, the failure of brittle materials, such as ceramics, is inseparable from the initiation and propagation of cracks. In this work, we will use the discrete element method to study the internal crack propagation mechanism of SiC ceramic materials under uniaxial compression. The effects of the pre-existing flaw, with inclination angles ranging from 0° to 75°, on the mechanical properties of the specimen were studied by the established DEM model of the SiC ceramic containing a single linear oblique flaw under uniaxial compression. At the same time, the evolution process of the force-chain field, displacement field and stress field around the pre-existing flaw in the process from loading to failure of the sample was also analysed.
2 . Methodology 2.1 Discrete element method Discrete element method is a numerical simulation method proposed by Cundall [15] in the early 1970s. It treats the material as a collection of discrete rigid particles which are connected by a bond. Potyondy and Cundall designed the bonded-particle model (BPM) [25] to simulate the mechanical behavior of bulk materials. In the BPM, it assumes that there exists a finite adhesive between two rigid particles, which can simultaneously transfer the force and moment between the particles. The assumed adhesive is represented by a parallel bond, and contacted particles are bonded together by means of parallel bonds to form an arbitrary shape assembly. The initiation and propagation of microcracks
in assembly can be seen as a breakage of the parallel bond between contacted particles. In this way, it can intuitively deal with the initiation and propagation process of crack, as shown in Fig. 1. In the DEM, the interaction between the contacted particles is a dynamic process, and the movement of particles at any time step meets Newton's second law. The movement has two ways of translation and rotation: Translation: Fi E g vi
Rotation: M i E g wi
m
I
'vi 't
'wi 't
(1)
(2)
where i (=1, 2) denotes the x and y coordinate directions, respectively; Fi is the out-of-balance force component of the particle; vi is the translational velocity; m is the mass of the particle; Mi is the out-of-balance moment due to the contacts; wi is the rotational velocity; I is the rotational inertia of the particle; βg is the global damping coefficient; Δt is the time step. Two arbitrary contacted particles are governed by the force-displacement rule, the contact force F is expressed as:
Fnn Fss
F
(3)
where Fn and Fs are the shear contact force and normal contact force respectively, the calculation formulas are as follows:
Fn 'F s
K nU n
(4)
K s 'U s
(5)
where Kn and Ks are normal and shear stiffness values at the contact, respectively; Un is the overlap of the two contacting particles; Fs is the increment of the shear contact force; ΔUs is the relative shear displacement (note that the shear contact force is initialised to zero when the two particles are in contact). The maximum tensile and shear stresses acting on the bond periphery are calculated (via beam theory) to be:
Vm a x
n M F R A I
(6)
s
Wm a x
Fi A
(7)
where A is the area of the cross section of a bond; I is the moment of inertia of the cross section of the bond; CF and M denote the total force and moment associated with the parallel bond respectively; CR is the radius of bond. The maximum tensile stress exceeds the tensile strength (σmax ≥ σbond) or the maximum shear stress exceeds the shear strength (τmax ≥ τbond), then the parallel bond breaks. And the relevant force, moment and stiffness are removed [27].
2.2 Modeling and calibration By PFC2D software developed by Itasca Company, two BPM models in which the particles were randomly distributed and closely arranged by adopting discrete element method was created in this paper. The dimensions of model are 4 × 8 mm and 2.4 × 0.5 mm, respectively, and the minimum radius of particles rlo was 1.5×10-2mm. Then, in order that the macroscopic mechanical properties of the model were consistent well with the actual SiC ceramic materials, as shown in Fig. 2, uniaxial compression test (the dimension of model is 4 × 8 mm), single edge notched beam test (the dimension is 2.4 × 0.5 mm, the span is 2 mm, the height of incision is 0.25 mm) and three-point bending test (the dimension is 2.4 × 0.5 mm, the span is 2 mm) were simulated respectively. By the trial and error method, the microscopic parameters of the discrete element model were constantly debugged and calibrated. The calibration results are shown in Table 1, it can be seen that the DEM simulation results are in good agreement with the experimental results. The discrete element model parameters of SiC ceramic are shown in Table 2. Table 1 Main mechanical properties of SiC in experimental tests and in DEM simulations Mechanical properties
Experimental results[21]
DEM simulation
Elastic Modulus (GPa)
420
422
0.14
0.14
Compressive strength (MPa)
2000
1935
Bending strength (MPa)
500-800
671
Fracture toughness (MPa/m1/2)
3.5
3.7
Poisson's ratio
Table 2 Parameters in the DEM model of SiC ceramic Variable name(unit)
Symbol
Value
Minimum radius of particles (m)
rlo
1.5e-5
Maximum-to-minimum radius ratio of the particle
n
1.2
Density of particle (kg/m3)
ρ
2600
Effective modulus of particle(Pa)
Eball
234e9
Normal-to-shear stiffness ratio of particle
kball
1.22
Effective modulus of parallel bond (Pa)
Ebond
234e9
Normal-to-shear stiffness ratio of parallel bond
kbond
1.22
Tensile strength of parallel bond (Pa)
σbond
9.3e8
Shear strength of parallel bond (Pa)
τbond
78e8
Friction coefficient of particles
μ
0.5
Deleted particles number of the pre-existing flaw
n
26
2.3 DEM model containing a single pre-existing flaw In order to produce a prefabricated central oblique flaw in the DEM model of SiC ceramic and ensure the width of resulting pre-existing flaw is as small as possible to meet the crack width distribution of the actual material at the same time, the method of deleting the particles was chosen. The diameter of the deleted particles represented the width of the flaw. And the half-length of the flaw can be obtained from:
a nu
d 2
n u rlo
(8)
where n is the deleted particles number along with flaw direction (as shown in Table 2), d is the minimum diameter of the particle, rlo is the minimum radius of the particle. In order to center the central oblique flaw with the center of the model, the length of the central oblique flaw can be determined by the deleted particles number along with flaw direction. The inclination angle θ of the center oblique flaw was represented by the included angle between the flaw direction and the positive direction of the x-axis. Moreover, θ can be arbitrarily specified as a controllable variable in the subsequent simulation. The model diagram is shown in Fig. 3.
2.4 Uniaxial compression test A uniaxial compression test of the discrete element model of SiC ceramics containing a single pre-existing flaw was established, as shown in Fig. 4. The pressure in the Y direction was applied to the specimen by setting a wall at the top and bottom of the discrete element model. The wall was represented by a horizontal line, whose stiffness was much larger than those of particles to prevent its deformation during the loading process. During the loading, both the top and bottom walls were moved to the center of the specimen with a constant v = 0.2 m/s by servo control. When the measured specimen stress is lower than 80% of the peak stress (defined as the compressive strength of the material), the loading was stopped and the specimen was broken finally. The inclination angle θ of the center oblique flaw in the specimen ranged from 0° to 75°. With 5° increment of inclination angle, those angles can be divided into 16 groups. So the fracture picture of the DEM model at different inclination angles was shown in Fig. 5. In those figures, the white straight line is the pre-existing center oblique flaw, while the red part indicates the microcracks generated in the loading process.
3. Results and Discussion
3.1 Analysis of Crack propagation During the loading process, the specimen went through four stages: initial compression, initiation of the first crack, generation of the secondary crack and specimen failure; for convenience, “stage 1”, “stage 2”, “stage 3” and “stage 4”, respectively, are used to represent these four stages. In Fig. 5 (stage 4), it was observed that the macroscopic crack angle was roughly between 60° and 70° when the specimen was broken at different flaw inclination angles. The failure modes of the specimen were basically the same. The break surface of the specimen that generated the macroscopic cracks produced many microcracks and resulted in many exfoliated blocks consisting of particles and exfoliated particles, as shown in Fig. 6. It is easy to infer that the largest stress is exerted on the specimen in the 60°~70° angle range. At the same time, there was a shear slip in this interval that finally resulted in a destabilizing rupture. Fig. 7 shows the initiation position of the first crack and secondary crack around the pre-existing flaw under compression. To eliminate the random error caused by the particle arrangement, we carried out six group numerical simulations on the specimen with the flaw inclination angles in the range of 0°~75° (interval of 5°). The enlarged images of the first crack and the secondary crack were chosen with the flaw angles at 0°, 25°, 45° and 65°, respectively, as shown in Figs. 8 and 9. It can be observed that the initiation position of the first crack was slightly different at the same flaw angle in Fig. 8. Additionally, the two growth tracks of the first crack were not symmetrical. With the increase of the flaw inclination angle θ, the distances d1 and d2 (the distance between the point of first crack initiation and flaw tip, respectively) were gradually reduced. At the same time, the two initiation angles of the first crack α1 and α2 increased (from approximately 90° to approximately 130°). These findings are in good agreement with the results of previous experiments [28] and simulations [16, 18]. Meanwhile, the initiation and propagation of the first crack were stable. A greater angle of the first crack makes it more difficult for the first crack to expand, and the propagation length of the crack was also gradually reduced. It has been widely confirmed in experimental studies of various brittle materials [30-33] that the initiation of the secondary crack at the flaw tips is normally subjected to shear force. Meanwhile, the secondary crack includes a coplanar secondary crack and an anti-wing crack, as shown in Fig. 7. The corresponding directions of the cracks were coplanar or nearly coplanar to the flaw and with an inclination similar to the wing cracks but in the opposite direction. Fig. 9 shows the enlarged images of the secondary crack when the flaw inclination angles are 0°, 25°, 45° and 65° (six distribution simulations were made for each flaw inclination angle.) It can be seen from the figure that the initiation position and direction of the secondary crack are in good agreement with the previous literature. The initiation and propagation of the secondary crack were unstable, and its crack surfaces were rough with some particles exfoliated by shearing.
3.2 Analysis of mechanical properties
The pre-existing flaw affected the mechanical properties of the specimen to a certain extent. Fig. 10 shows the DEM simulation results of the mechanical properties of SiC ceramics with a flaw for which the inclination angle ranged from 0° to 75° (the results were the average of the results of six different simulations, and the error bars represents the relevant standard deviation). The horizontal and vertical axes represent the flaw inclination angle and the ratio of the crack initiation stress to the compressive strength of the specimen, respectively. By comparing to the experimental results of for gypsum [28] and the DEM simulation results for a rock-like material [17], it can be concluded that the corresponding trend was roughly the same. At the same time, the ratios in the rock-like material and the SiC ceramics in the DEM simulation results were smaller than that obtained experimentally for the Gypsum material. One plausible explanation is that this was due to experimental limitations. In the physical tests, the detection of the first crack was based on the first distinct cracking sound or the minute crack opening event observed on the specimen surface. As suggested in other experimental studies employing the AE and CT techniques, cracking events prior to those observed on the surface could have occurred inside the specimens. A proper identification of the latter could lead to the lower ratio [17]. Fig. 11 shows the trends in Young's modulus and Poisson's ratio of the specimen when the flaw inclination angle ranged from 0° to 75°. It can be seen that the presence of the flaw in the specimen reduces the effective Young's modulus (420 GPa with no flaw) of the specimen and increases the Poisson's ratio (0.14 with no flaw) obtained in the uniaxial compression test. With the increase of the flaw inclination angle, the effective Young's modulus changed from being relatively stable to increasing, but the obtained value was always smaller than the original value. Meanwhile, Poisson's ratio was first relatively stable and then decreased, but it always remained higher than its original value.
3.3 Analysis of force-chain Fig. 12 shows the distribution of the force-chain field around the flaw at four different stages for the flaw inclination angles of 0°, 25°, 45° and 65° (stages 1, 2, 3 and 4 are the processes of initial compression, initiation of the first crack, generation of a secondary crack and specimen failure, respectively). In the figure, the particles are hidden, and the green and the black lines represent the tensile and the compressive force-chains of the parallel bond, respectively. Additionally, the thickness of the force-chain lines represents the magnitude of the contact resultant force on the parallel bond, while its direction indicated the direction of the contact resultant force. Finally, the red represented the microcracks generated during the loading process. Because in the DEM simulation the generation of microcracks is characterized by the breaking of the parallel bonds between the particles, the magnitude and direction of the force on the parallel bond around the flaw in the loading process (the distribution of the force-chain field) can, to a certain extent, reflect the force distribution and evolution around the flaw. Fig. 13 shows the evolution process of the strong force-chain field around the flaw for the inclination angles of 0°,
25°, 45° and 65° at different loading stages. Combining Figs. 12 and 13, we can observe that at the initial compression stage (stage 1), the concentration regions of the tensile force-chain (the green region A) were generated in the top and bottom of the pre-existing flaw at different inclination angles. Meanwhile, the concentration regions of the compressive force-chain (the black region B) were generated in the left and right edges of the pre-existing flaw. This finding indicated that a stress concentration was generated in regions A and B. However, the position and size of the concentration regions of the tensile force-chain (the green region A) changed for different flaw inclination angles. For example, for the flaw inclination angle of 0°, the A regions were just above or below the flaw, and the A region area was the largest. With the increase of the flaw inclination angle, the regions were gradually reduced and transferred to the flaw tip. By comparing the difference of the first crack initiation position in section 3.1 above, we can draw the conclusion that the first crack initiation position occurred in the concentration regions of the tensile force-chain (the green region A). With the decrease of the flaw inclination angle, both the concentration regions of the tensile force-chain and the distribution interval of the first crack initiation position become larger. Thus, the decrease of the flaw inclination angle in the experiments reported in the literature [29] increased the distribution error of the distance d1 and d2 between the first crack initiation position and the flaw tip instead. With the constant loading, the first crack in stage 2 was generated in the concentration regions of the tensile force-chain (the green region A) and then was consistently extended. Concentrated tensile stress was released during the first crack initiation and extension. Meanwhile, the region A transferred to the first crack tip and gradually shrank. When the secondary crack occurred in stage 3, the first crack stopped extending, and the tensile force-chain concentration region (the green region A) disappeared. However, this process did not have a significant effect on the compressive force-chain concentration region (the black region B) at both ends of the flaw. In addition, around the flaw, a new region called the “force-chain dissipation region” (the blue region C) was generated by the initiation of the first crack and changed accordingly with the extension of the first crack. It can be seen from Fig. 13 that the increasing flaw inclination angle made it difficult for the first crack to extend. As a result, the length of the first crack length decreased and the area of region C decreased as well. To a certain extent, the “force-chain dissipation region” (the blue region C) ensured the steady extension of the first crack as well as the absence of a new crack branch in the process of first crack extension. A quantitative analysis of the force-chain in the specimen was conducted in the following figures. Fig. 14 shows the distribution between the normal force magnitude of the parallel bond and the force-chain angle at different loading stages when the flaw inclination angles were 0°, 25°, 45° and 65°. The positive value of the normal force denoted compression while the negative value denoted tension. The force-chain angle was the included angle between the normal unit vector direction of the parallel bond and the positive direction of the x-axis. In Figs. 14(a-d), in the initial stage, the normal force of the parallel bond was initialized to zero. As the loading was carried out, the normal force
began to increase, and the top and bottom of the corresponding curve also began to extend in the 90° and 270° directions, indicating the maximums of the normal compression force on the parallel bond around the angle intervals at that time. However, the left and right ends of the curve were retracted in the opposite direction towards 180° and 0°, indicating that in the angle intervals, the normal tension force on the parallel bond was maximized. Finally, the final curve gradually became a “peanut” rose diagram. When the specimen failed in stage 4, the value of the normal force on the parallel bond showed a drop at the different flaw inclination angles, and the curve moved back by a certain amount. Meanwhile, it also found that the curve was sunken slightly in the 90° and 270° direction, and the normal force on the parallel bond decreased slightly in the specimen. It is possible that the propagation and coalescence of the crack in the specimen weakens the load-bearing capacity of the specimen at the angle intervals. Fig. 15 shows the distribution between the shear force magnitude of the parallel bond and the force-chain angles at the different loading stages in the specimen for the flaw inclination angles of 0°, 25°, 45° and 65°. As seen from Figs. 15(a-d), the shear force on the parallel bond is initialized to zero. As the load increased, the shear force increased accordingly. The curves in the angle intervals in the 50°-55°, 140°-145°, 230°-235° and 320°-325° intervals continue to extend, and the shear force reached maximum values in the four angle intervals. As a result, the final curve gradually became the "butterfly wings" rose diagram. When the specimen failed in stage 4, the shear force magnitudes of the parallel bond dropped at the different flaw inclination angles, and the curve moved back by a certain amount back as well. The comparison between the above study and the angle interval range of the generated macrocracks in the previous section 3.1 confirmed the internal origin of the macroscopic cracks generated in the 60° ~ 70° angle interval when the specimen failed.
3.4 Analysis of displacement During the process from initial compression to rupture failure of the discrete element specimen, the displacement of the particles can be indicated by hiding the particles and only showing the displacement tendency using arrows. To some extent, this facilitated the observation of the process of microcrack initiation, propagation and coalescence at different loading stages. The displacement tendency arrows show a corresponding evolution at different loading stages. As a result, we divided the displacement field into four different types, denoted type I, type II, type III, and type IV for convenience. As shown in Fig. 16, type I represents a purely tensile displacement field, type II represents both the tensile and shear displacement field, type III represents a displacement field that was neither tensile nor shear, and type IV represents both the tensile and shear displacement fields but with the tensile field much larger than the shear field. Figs. 17 (a-d) show the evolution process of the displacement tendency arrows for the flaw inclination angles of 0°, 25°, 45° and 65°. It can be seen that the arrows around the first crack were open to both sides of the first crack regardless of the change of the flaw inclination angles at stage 2. At that time, the displacement field of the first crack
was type Ϩ, indicating that the first crack was generated by tension. As the loading progressed, the secondary cracks in stage 3 were generated accordingly. The displacement tendency arrows around the crack initiation position gradually moved from their original upward or downward positions towards the direction far away from both sides of flaw tips after the arrows approached the flaw tips. The displacement field of the secondary cracks was almost identical to type II, which meant that the secondary crack was generated by both tension and shear. At the same time, we also observed that the displacement tendency arrows around the first crack began to evolve from type Ϩ to type Ϫ because of secondary crack generation and some of first cracks gradually stopped extending. During the constant loading process, some of the secondary cracks gradually developed into macroscopic cracks. When the specimen was completely broken in the stage 4, the displacement field showed more changes. The changes were as follows: the displacement tendency arrows around the first crack almost turned into type III because the first crack no longer expanded while the arrows around the secondary crack shifted from the original type ϩ to type Ϫor type ϫ. The displacement field where the secondary crack evolved into the macroscopic crack changed from the original type II to type ϫ, and the displacement field where the secondary crack stopped the initiation transformed into type III.
3.5 Analysis of stress The study of the evolution of the stress field around the flaw at different loading stages is beneficial for further understanding of the entire process of crack initiation, propagation and convergence. By setting up a series of measurement circles around the flaw, the measurement of the stress distribution around the flaw was carried out in this work and then a stress distribution cloud picture was drawn. The measured stress magnitude was the average of the stress carried by the particles contained in the circle. The number of particles was generally 4-6. The distribution of the measurement circle is shown in Fig. 18. Here, the stress distributions in stages 1, 2 and 3 were selected and measured for the crack inclination angles of 0°, 25°, 45° and 65°. Then, the stress distribution cloud picture of the three stress components, namely, σxx, σxy and σyy, were generated accordingly. As shown in Figs. 19-21, the coordinate system in the stress cloud picture adopted a local coordinate system of the flaw. The centre of the flaw was the origin of the coordinate system, and the x- and y-coordinate values divided by the half-length of the preset crack were the dimensionless numbers. A negative stress value denoted compression and a positive value corresponded to tension, in the units of MPa. Figs. 19(a-d) showed the σxx stress component distribution around the flaw at different stage for the flaw inclination angles of 0°, 25°, 45° and 65°. For the flaw inclination angle of 0°, we observed that the compressive stress concentration regions (blue regions) appeared at the flaw tips. Meanwhile, a large range of tensile stress concentration regions (red regions) appeared immediately above and below the flaw. With the increase of the flaw inclination angle, the compressive stress concentration regions at the flaw tips moved correspondingly to the centre of the flaw.
Meanwhile, the tensile stress concentration regions above and below the flaw also moved towards the flaw tips. In the process of changing the inclination angle of the flaw from 0° to 65°, we found that both the compressive stress concentration regions and the tensile stress concentration regions around the flaw exchanged their positions. With the loading, the first crack was generated and extended in the stage 2. It can be observed in the figure that the first crack initiated in the tensile stress concentration regions (red regions) and released the tensile stress in this region. With the extension of the first crack, the original tensile stress concentration regions transferred to the first crack tips, while the compressive stress concentration regions (blue regions) were not affected. Then, the secondary crack initiated and extended near the flaw tips in stage 3, accelerating the dissipation of the tensile stress concentration regions. However, there was no effect on the tensile stress concentration regions. Fig. 20 shows the σxy stress component distribution around the flaw with different inclination angles at different stages. The transfer phenomenon of the stress concentration regions was also observed in stage 1 when the specimen was loaded at the same time. In the process of changing the inclination angle of the flaw from 0° to 65°, the positive stress concentration regions were initially located in the left bottom and right top of the flaw rotated clockwise 180° around the flaw tip and then transferred to the left top and right bottom of the flaw. The negative stress concentration region was initially located in the left top and right bottom of the flaw and rotated clockwise around the flaw and transferred to the right top and left bottom of the flaw. It is important to note that the initiation and expansion of the first and secondary cracks in stages 2 and 3 did not affect the evolution of the σxy stress component. Fig. 21 shows the distribution of the stress component σyy around the corresponding flaw. The stress component σyy in the stress cloud picture was negative due to the compression loading of the specimen. It was obvious that the compressive stress concentration regions were generated at the flaw tips (blue regions). Meanwhile, it appears shielding regions of compressive stress (brown regions) were formed at the top and bottom of the flaw. The shielding effect was stronger at smaller distances from the flaw surface. As a result, the compressive stress became weaker.
3.6 Discussion and outlook As shown in Figs. 10 and 11, it is important to note that when the flaw inclination angle was in the range of 30-40º, all curves of the ratio of crack initiation stress to compressive strength, Young's modulus and Poisson's ratio changed. The possible reasons for the change are as follows: (1) for the length of first crack, when the flaw inclination angle was in the range of 0-30º, the first crack was easier to extend. This characteristic resulted in the greater length of the first crack when the specimen failed under loading (as shown in Figs. 8 and 9). However, for the flaw inclination angles greater than 30º, it became more difficult for the first crack to extend. Therefore, the length of the first crack for the flaw inclination angles greater than 30º was relatively shorter. The difference in the length of the first crack resulted in the difference in the energy consumption at the flaw tips. Accordingly, the mechanical properties of the specimen were
different. (2) For the size and dissipation of the concentration regions of the tensile force-chain and tensile stress, when the flaw inclination angle was in the range of 0-30º, the concentration regions of the tensile force-chain (green regions in Figs. 12 and 13) and the tensile stress (red regions in Fig. 19) near the flaw tips were larger. The tensile stress concentration phenomenon at the flaw tips was more distinct than for the cases where the flaw inclination angle was greater than 30º. As a result, when the first crack was generated and extended, the “force-chain dissipation region” was larger, and the tensile stress dissipation effect was more obvious. In this work, a model of a SiC ceramic containing a single pre-existing flaw was established based on the discrete element method. However, the distribution of the microcrack is random and their direction cannot be controlled in the process of preparation, processing and service of the actual material. For engineering ceramics, the microcracks caused by machining damage and intrinsic defects can be detected and measured by non-destructive testing techniques (penetrant testing, ultrasonic testing, laser ultrasonication and infrared testing). The different flaw inclination angles will result in the differences of the mechanical properties of the specimens. To some extent, we can predict and evaluate the mechanical properties of the specimen using the DEM modelling based on the reasonable statistics for the distribution of the microcrack with different inclination angles in the specimen. Thus, this work can increase the reliability and service life of the specimen under different service conditions. At the same time, ceramic materials show a strong sensitivity to crack defects and a high fragility. An unpredictable brittle fracture in service is mostly caused by the initiation, propagation and coalescence of interval cracks. Ceramic toughening techniques can effectively improve the fracture toughness of ceramics as well as their reliability and service life. To enhance the toughness of ceramic materials, it is crucial to improve their ability to resist crack propagation and to reduce the stress concentration effect of crack tips. Studies of the relationship between the flaw inclination angle and the location and angle orientation of the toughening material in specimen can further reveal the inherent mechanism of ceramic toughening and explore the effect of the location and arrangement of toughening materials around the flaw on the toughening effect of the specimen under the conditions of different flaw inclination angles. Additionally, such studies can provide guidance for the efficient toughening of ceramic materials.
4
Conclusion (1) The flaw reduces the effective Young's modulus of the specimen and increases the Poisson's ratio under the
uniaxial compression test. With the increase of the flaw inclination angle, the effective Young's modulus changes from a being relatively stable to increasing, but its value is always smaller than the original value. Meanwhile, Poisson's ratio is first relatively stable and then decreases, and this value is always greater than its original value. (2) In the process of changing the inclination angle of the flaw from 0° to 65°, the concentration regions of the tensile force-chain above and below the flaw are gradually reduced and transferred to the flaw tip. The first crack
generated in the concentration regions of the tensile force-chain makes the region transfer to the first crack tip and decrease gradually. (3) During the loading process, the distribution curve between the normal force magnitude of the parallel bond and the force-chain angle gradually transforms into the "peanut" rose diagram. The shear force curve with angles gradually becomes similar to the "butterfly wings" rose diagram. When the specimen fails, the force magnitude of the parallel bond shows a drop. The study of the normal and shear forces of the parallel bond with the angle distribution in the specimen confirms the intrinsic origin of the macroscopic crack generation in the 60°-70° angle interval. (4) The evolution of the displacement field around the flaw in the process from the first crack initiation to the specimen failure is analysed. The initiation, propagation and coalescence of the crack are summarized by four types of displacement tendency arrows (type I, type II, type III and type IV). It is confirmed that the initiation of the first crack is affected by the tensile force while the secondary crack is affected by both the tensile and shear forces. (5) The first crack initiated in the tensile stress concentration regions. With the extension of the first crack, the tensile stress in the region is released. The tensile stress concentration regions transfer to the first crack tip and decrease. However, the initiation and propagation of the first and secondary cracks have little effect on the compressive stress concentration regions at the flaw tips. At the same time, examination of stress cloud picture analysis of the σxx and σxy stress components around the flaw under initial compression shows that the positive and negative stress concentration regions are transferred clockwise along the flaw in the process of changing the flaw inclination angles from 0° to 65°. In the stress cloud picture analysis of the σyy stress component, the compressive stress concentration is generated at the flaw tips. It appeared to form shielding regions of compressive stress at the top and bottom of the flaw. Stronger shielding effects were observed closer to the flaw surface. As a result, the compressive stress became weaker.
Acknowledgement This work was supported by the National Natural Science Foundation of China (Grant No. 51605409, 11772135, 51475404 and 11602212)
References [1] B. Lawn, Fracture of brittle solids, Cambridge university press, Cambridge, 1993. [2] H. Abe, M. Naito, T. Hotta, Flaw size distribution in high-quality alumina, J. Am. Ceram. Soc. 86 (2010) 1019-1021. [3] W.H. Tuan, J.C. Kuo, Effect of abrasive grinding on the strength and reliability of alumina, J. Eur. Ceram. Soc. 18 (1998) 799-806. [4] K. Ikeda, H. Igaki, Effect of surface flaw size on fracture strength of alumina ceramics, J. Am. Ceram. Soc. 70
(1987) C-29-C-30. [5] Z.T. Bieniawski, Mechanism of brittle fracture of rock: parts 1, Int. J. Rock. Mech. Min. 4 (1967) 395-404. [6] L.N. Germanovich, R.L. Salganik, A.V. Dyskin, Mechanisms of brittle fracture of rock with pre-existing cracks in compression, Pure. Appl. Geophys. 143 (1994) 117-149. [7] A. Bobet, H.H. Einstein, Fracture coalescence in rock-type materials under uniaxial and biaxial compression, Int. J. Rock. Mech. Min. 35 (1998) 863-888. [8] P. Cao, T. Liu, C. Pu, Crack propagation and coalescence of brittle rock-like specimens with pre-existing cracks in compression, Eng. Geol. 187 (2015) 113-121. [9] C.H. Park, A. Bobet, Crack coalescence in specimens with open and closed flaws: a comparison, Int. J. Rock. Mech. Min. 46 (2009) 819-829. [10] B. Shen, The mechanism of fracture coalescence in compression-experimental study and numerical simulation, Eng. Fract. Mech. 51 (1995) 73-85. [11] L.N.Y. Wong, H.H. Einstein, Crack coalescence in molded gypsum and carrara marble: part 1-macroscopic observations and interpretation, Rock. Mech. Rock. Eng. 42 (2008) 475-511. [12] L.N.Y. Wong, H.H. Einstein, Crack coalescence in molded gypsum and carrara marble: part 2-microscopic observations and interpretation, Rock. Mech. Rock. Eng. 42 (2009) 513-545. [13] T.J. Ahrens, Rock physics & phase relations: A handbook of physical constants, American Geophysical Union, Washington DC, 1995. [14] X.T. Feng, S. Chen, H. Zhou, Real-time computerized tomography (CT) experiments on sandstone damage evolution during triaxial compression with chemical corrosion, Int. J. Rock. Mech. Min. 41 (2004) 181-192. [15] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Géotechnique. 29 (1979) 47-65. [16] X.P. Zhang, L.N.Y. Wong, Crack initiation, propagation and coalescence in rock-like material containing two flaws: a numerical study based on bonded-particle model approach, Rock. Mech. Rock. Eng. 46 (2013) 1001-1021. [17] X.P. Zhang, Cracking processes in rock-like material containing a single flaw under uniaxial compression: a numerical study based on parallel bonded-particle model approach, Rock. Mech. Rock. Eng, 45 (2012) 711-737. [18] L.A.M. Camones, E. do Amaral Vargas, R.P.D. Figueiredo, Application of the discrete element method for modeling of rock crack propagation and coalescence in the step-path failure mechanism, Eng. Geol. 153 (2013) 80-94. [19] S. Hentz, F.V. Donzé, L. Daudeville, Discrete element modeling of concrete submitted to dynamic loading at high strain rates, Comput. Struct. 82 (2004) 2509-2524. [20] M.Y. Zang, Z. Lei, S.F. Wang, Investigation of impact fracture behavior of automobile laminated glass by 3D discrete element method, Comput. Mech. 41 (2007) 73-83.
[21] Y.Q. Tan, D.M. Yang, Y. Sheng, Discrete element method (DEM) modeling of fracture and damage in machining process of polycrystalline SiC, J. Eur. Ceram. Soc. 29 (2009) 1029-1037. [22] Y.Q. Tan, D.M. Yang, Y. Sheng, Study of polycrystalline Al2O3 machining cracks using discrete element method, Int. J. Mach. Tool. Manu. 48 (2008) 975-982. [23] Y.Q. Tan, S.Q. Jiang, D.M. Yang, Scratching of Al2O3, under pre-stressing, J. Mater. Process. Tech. 211 (2011) 1217-1223. [24] S.Q. Jiang, X. Li, Y.Q. Tan, Discrete element simulation of SiC ceramic with pre-existing random flaws under uniaxial compression, Ceram. Int. 43 (2017) 13717-13728. [25] D.O. Potyondy, P.A. Cundall, A bonded-particle model for rock, Int. J. Rock. Mech. Min, 41 (2004) 1329-1364. [26] N. Cho, C.D. Martin, D.C. Sego, A clumped particle model for rock, Int. J. Rock. Mech. Min. 44 (2007) 997-1010. [27] P.A. Cundall, PFC2D user's manual (Version 3.1), Minnesota: Itasca Consulting Group Inc, Minnesota, 2004. [28] L.N.Y. Wong, H. Einstein, Fracturing behavior of prismatic specimens containing single flaws, Golden Rocks 2006, The 41st US Symposium on Rock Mechanics (USRMS), American Rock Mechanics Association, 2006. [29] A. Bobet, The initiation of secondary cracks in compression, Eng. Fract. Mech. 66 (2000) 187-219. [30] E.Z. Lajtai, Brittle fracture in compression, Int. J. Fracture. 10 (1974) 525-536. [31] A.R. Ingraffea, F.E. Heuze, Finite element models for rock fracture mechanics, Int. J. Numer. Anal. Met. 4 (1980) 25-43. [32] H. Jiefan, C. Ganglin, Z. Yonghong, An experimental study of the strain field development prior to failure of a marble plate under compression, Tectonophysics. 175 (1990) 269-284. [33] Y.P. Li, L.Z. Chen, Y.H. Wang, Experimental research on pre-cracked marble under compression, Int. J. Solids. Struct. 42 (2005) 2505-2516.
List of figure captions: Fig. 1. BPM model diagram in Ref.[26]. Fig. 2. Calibration tests of mechanical properties, (a) uniaxial compression test, (b) single edge notched beam test, and (c) three-point bending test. Fig. 3. (a) Model with a single flaw, (b) DEM model containing a single pre-existing flaw. Fig. 4. DEM model in uniaxial compression test. Fig. 5. Rupture picture of specimens at different flaw inclination angles (0°~75°). Fig. 6. (a) Macroscopic crack in specimen with 45° pre-existing flaw, (b) enlarged image of yellow circle “A” and (c) enlarged image of yellow circle “B”. Fig. 7. Initiation position diagram of the first crack and secondary crack (θ is the inclination angle of pre-existing flaw. α1 and α2 are the crack initiation angles. d1 and d2 are the distances between the point of first crack initiation and flaw tip)[28,29]. Fig. 8. Enlarged image of the first crack under different flaw inclination angles, (a) 0°, (b)25°, (c)45° and (d)65°. Fig. 9. Enlarged image of the secondary crack under different flaw inclination angle, (a) 0°, (b) 25°, (c) 45° and (d) 65°. Fig. 10. Comparison of the ratio of the crack initiation stress and compressive strength under different flaw inclination angles. Fig. 11. Young’s modulus and Poisson’s ratio under the different flaw inclination angles. Fig. 12. Distribution of the force-chain field around the flaw at four stages under the different flaw inclination angles (a) 0°, (b) 25°, (c) 45° and (d) 65° (green lines indicate the tensile force-chain, black lines indicate the compressive force-chain, and red lines indicate the ruptured parallel bonds). Fig. 13. Distribution evolution of the strong force-chain field around the flaw at different loading stages when the flaw inclination angles are (a) 0°, (b) 25°, (c) 45° and (d) 65° (green regions A indicate the concentration regions of the strong tensile force-chain, black regions B indicate the concentration regions of the strong compressive force-chain, and blue regions C indicate the “force-chain dissipation region”). Fig. 14. Contrast analysis of the normal force of parallel bond vs the angle under different flaw inclination angles. Fig. 15. Contrast analysis of the shear force of parallel bond vs the angle under different flaw inclination angles. Fig. 16. Four types of displacement tendency arrows at different loading stages. Fig. 17. Evolution process of displacement tendency arrows at different loading stages with the flaw inclination angles (a) 0°, (b) 25°, (c) 45°, and (d) 65°. Fig. 18. Measurement circle distribution diagram around the pre-existing flaw. Fig. 19. Stress component σxx around the flaw at different loading stages with different inclination angles: (a) 0°, (b) 25°, (c) 45°, and (d) 65°.
Fig. 20. Stress component σxy around the flaw at different loading stages with different inclination angles: (a) 0°, (b) 25°, (c) 45°, and (d) 65°. Fig. 21. Stress component σyy around the flaw at different loading stages with different inclination angles: (a) 0°, (b) 25°, (c) 45°, and (d) 65°.
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH
)LJXUH