Discrete element simulation of SiC ceramic with pre-existing random flaws under uniaxial compression

Author’s Accepted Manuscript Discrete element simulation of SiC ceramic with pre-existing random flaws under uniaxial compression Shengqiang Jiang, Xu Li, Yuanqiang Tan, Haohan Liu, Zhiqiang Xu, Rui Chen www.elsevier.com/locate/ceri

PII: DOI: Reference:

S0272-8842(17)31516-X http://dx.doi.org/10.1016/j.ceramint.2017.07.084 CERI15804

To appear in: Ceramics International Received date: 19 June 2017 Revised date: 8 July 2017 Accepted date: 11 July 2017 Cite this article as: Shengqiang Jiang, Xu Li, Yuanqiang Tan, Haohan Liu, Zhiqiang Xu and Rui Chen, Discrete element simulation of SiC ceramic with preexisting random flaws under uniaxial compression, Ceramics International, http://dx.doi.org/10.1016/j.ceramint.2017.07.084 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 with pre-existing random flaws under uniaxial compression Shengqiang Jianga*, Xu Lia, Yuanqiang Tanb, Haohan Liua, Zhiqiang Xua, 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. [email protected]

Abstract Because of the outstanding advantages in studying crack propagation and coalescence of brittle solids, the discrete element method (DEM) has been widely used to study the crack propagation in the processing of engineering ceramics. To explore the effects of random defects on the failure mode and mechanical properties of an SiC ceramic material, a DEM model of the SiC ceramic material was established using the Mori–Tanaka method wherein the crack density was determined using the number and length of pre-existing flaws. The established DEM model with the pre-existing random flaws of the SiC ceramic was then subjected to a uniaxial compression test. To eliminate the randomness of the established model as much as possible, more DEM tests were performed with five different randomness flaw distributions. With the increase in the crack density, the compressive strength and crack initiation stress of the material significantly declined. Based on both the fracture specimen and stress–strain curve, the failure mode of the specimen could be divided into two stages: brittle fracture and plastic failure. In the interval with crack densities lower than 0.028, the specimen ruptured along angles in the range of 55°–60°, and the failure mode was brittle fracture. Moreover, the effective Young’s modulus of the specimen obtained using the DEM was in good agreement with that obtained using the Mori–Tanaka method curve. However, when the crack density was higher than 0.028, which is regarded as a high crack-density interval, the cracks propagated and coalesced in a direction parallel to that of the maximum principal stress. In addition, the failure mode was largely axial splitting. The effective Young’s modulus obtained using the DEM was lower than that obtained using the Mori–Tanaka method, wherein only the weak interaction between the cracks was considered. Keywords: Pre-existing random flaws；Crack-density；SiC；Mechanical property；Discrete element method 1

Introduction Ceramic materials with excellent mechanical, thermal, and chemical properties are widely used in machinery, aerospace,

microelectronics, and other important fields [1]. However, in ceramics, which is as a typical brittle material, various random defects are induced in the preparation, processing, and service of materials, thereby resulting in uncontrollable effects with regard to material strength and reliability. Microcracks are a common type of internal defects in brittle materials. The study on its mesoscopic damage and rupture mechanism has been a hot topic, which is significant in engineering applications for material performance testing, evaluation, and research [2, 3]. The studies on mesoscopic damage and rupture mechanism in microcracks are very complicated. The damage induced by crack defects is an essential mechanism in brittle solids subjected to a compressive stress. Ashby et al. [4] analysed the growth of a crack from an initially inclined flaw and examined the conditions under which an array of such cracks interacts. Cao et al. [5] investigated the crack propagation and coalescence by loading rock-like specimens with two and three pre-existing flaws, which were induced by pulling out the embedded metal inserts in the pre-cured period. Kushch et al. [5, 6] studied the effects of angle distribution and crack density in cracked solids on the effective rigidity. Several studies have been conducted on research methods for materials with microcrack defects, and methods of calculating the effective elastic modulus of the materials with microcrack defects have received increasingly more attention [8-11]. Mori and Tanaka proposed the Mori–Tanaka method by considering the weak interaction between the microcracks [12]. Benveniste [13] studied the application of the Mori–Tanaka method in calculating the effective modulus of microcracks. In the above

·2·

studies, the brittle solids with crack defects were largely analysed via experiments and theoretical analysis. However, it is difficult to manufacture specimens containing defects, particularly the ones with random defects, because of the limited processing capabilities. Moreover, the theoretical analysis based on the fracture mechanics cannot appropriately deal with the random distribution and various shapes of defects. Therefore, it is necessary to explore a new and reliable numerical method to investigate the effects of different crack densities on the mechanical properties of the material containing random defects. In this study, a discrete element method (DEM) model of an SiC ceramic is established. Based on the Mori–Tanaka method, the crack-density function, wherein the number and length of pre-existing flaws are variables, is used to develop the DEM model of the SiC ceramic with pre-existing random flaws. The failure mode and mechanical properties of the SiC ceramic specimen with different crack densitie an SiC ceramic s are analysed under uniaxial compression. Moreover, the effective Young’s modulus obtained using the DEM is compared with the theoretical curve obtained using the Mori– Tanaka method. 2

Discrete element method Cundall [7] proposed the DEM, which is a numerical simulation method, in the early 1970s to study the rock mechanics

and was later applied to study the fracture, failure, and processing of brittle materials such as rocks [15,16], concrete [17], laminated glass [18,19], and ceramics [20-23]. In the DEM, the continuous bulk material is considered a collection of discrete particles connected by bonds. The microcrack initiation and propagation are regarded as the breakage of bonds between the particles; hence, it has outstanding advantages in terms of the initiation and propagation of microcracks, as demonstrated in Fig. 1. In the DEM, the interaction between the contacting particles is treated as a dynamic process, and the stress and deformation of the entire particle assembly are obtained from tracing the movement and interactions of each individual particle. Each contact between the particles can be described using springs, friction resistance, and damp absorber, as shown in Fig. 2. The solution of a DEM model is obtained by applying Newton’s second law and the force –displacement law periodically to form one cycle of computation. Newton’s second law is used to calculate the acceleration of the particle resulting from the two contacting elements, whereas the force–displacement law is used to update the contact forces resulting from the overlap of the two contacting elements. The two laws are applied repeatedly to obtain the entire calculation cycle of the DEM. The motion of the particles is governed by Newton’s second law as follows: Translation motion

Rotation motion

vi t

(1)

w M i   g w i I i t

(2)

Fi   g vi  m

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. For the force–displacement law, the resulting contact-force vector at the contact between any two particles F can be resolved into normal and shear components Fn and Fs as follows.

F  Fnn  Fss The magnitudes of the normal and shear contact forces are calculated as follows:

(3)

n

Fn  Kn Un

(4)

F s   K s U s

(5)

s

where K and K 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 stiffness Ki (i=n, s) at the contact is obtained from the following equation.

Ki 

kiAkiB kiA  kiB

(6)

where kiA and kiB (i = n, s) are the normal and shear stiffness values of the two contacting particles, respectively, which can be calculated by knowing the ball effective modulus Eball and ball normal-to-shear stiffness ratio kball. The contact is checked for slip conditions by calculating the maximum allowable shear-contact force as follows.

Fs   Fn

(7)

where μ is the friction coefficient at the contact. The force–displacement law is used to analyse the forces in parallel bonds. Fig. 3 shows a schematic of the parallel-bond contacts[24]. Here, it is assumed that there is a finite size of adhesive represented by the parallel bond between two rigid particles, which can simultaneously transfer the forces and moments between the particles. The total force and moment associated with the parallel bond are denoted by Fi and M, respectively. The force and moment vectors can be resolved into normal and shear components with respect to the contact plane as follows: n

s

Fi  Fi  Fi n

(8) s

Mi  Mi Mi

(9)

where Fin andFis andMin andMis denote the normal and shear force and moment component vectors, respectively, on the contact plane. When the bond is formed,Fi andMi are initialised to zero. Each subsequent relative displacement and rotation increment at the contact results in an increment in the elastic force and moment, which is added to the current values. The elastic-force increments occurring over a time step of Δt are calculated as follows. n

n

 F i  (k AU n )ni s

s

 F i  k AU is with

(10) (11)

U i  vi t

(12)

The elastic-moment increment is calculated as follows: n

 M  k I  with

(13)

  (w   w  )t B

A

(14)

where kn and ks denote the normal and shear stiffness values of the parallel bond, which can be calculated by knowing the bond effective modulus Ebond and bond normal-to-shear stiffness ratio kbond; A is the area of the cross section of a bond; I is the moment of inertia of the cross section of the bond about an axis through the contact point and in the direction of Δθ. The maximum tensile and shear stresses acting on the bond periphery are calculated (via beam theory) as follows.

·4·

 max

n M F   R A I

(15)

s

 max 

Fi

(16)

A

If the maximum tensile stress exceeds the normal strength (σmax ≥ σbond) or the maximum shear stress exceeds the shear strength (τmax ≥ τbond), the parallel bond breaks. The relevant force, moment, and stiffness are then removed. When the parallel bond breaks continually along a certain direction, microcracks generate and propagate along that direction. Fig. 4 shows the constitutive behaviour of the bond [24, 25]. 3

DEM model of SiC ceramic To establish a DEM model similar to the real SiC ceramic material, dense packing particles were arranged closely and

randomly, generated by using the PFC 2D software developed by Itasca company [24]. The length w, height l, and minimum-particle radius rlo of the model were 4 mm, 8 mm, 0.015 mm, respectively. The bonded particle model (BPM) [26] is a particle assembly with arbitrary shape wherein the particles are connected together by a parallel bond and will be used to simulate the mechanical behaviours of the solid materials. By debugging and resetting the microscopic parameters of the particles and parallel bonds in the model repeatedly, the macroscopic mechanical properties of the SiC ceramic in the DEM model are found to be consistent with the actual material. In the process, the parameters pertaining to the mechanical properties of the established model were calibrated by conducting uniaxial compression test, single-edge notched beam test, and three-point bending test, as shown in Fig. 5. Table 1 lists the DEM-simulation results of the mechanical properties, such as compressive strength, Young’s modulus, and Poisson’s ratio in the uniaxial compression test, bending strength in the three-point bending test, and fracture toughness in the single-edge notched beam test. They are in good agreement with the experimental results. Table 2 lists the parameters of the SiC ceramic in the DEM model. Table 1 Main mechanical properties of SiC in experimental test and in DEM simulations Mechanical properties Experimental results [21] DEM simulation Elastic Modulus E(GPa)

420

422

Poisson's ratio

0.14

0.14

Compressive strength (MPa)

2000

1935

Bending strength (MPa)

500-800

671

Fracture toughness(MPa/m1/2)

3.5

3.7

Table 2 Parameters in the SiC-ceramic discrete element model Variable name(unit) Symbol

Value

Length of model(m)

w

4e-3

Height of model(m)

l

8e-3

Minimum radius of ball (particles)(m)

rlo

1.5e-5

ball maximum-to-minimum radius ratio

n

1.2

Density of model(kg/m3)

ρ

2600

Ball effective modulus(Pa)

Eball

234e9

Ball normal-to-shear stiffness ratio

kball

1.22

Bond effective modulus(Pa)

Ebond

234e9

4

Bond normal-to-shear stiffness ratio

kbond

1.22

Bond tensile strength(Pa)

σbond

9.3e8

Bond shear strength(Pa)

τbond

78e8

Friction coefficient

μ

0.5

DEM model of SiC ceramic with pre-existing random flaws

4.1 Pre-existing random flaws The Mori–Tanaka method is a method for calculating the effective modulus of a solid with microcracks. The main principle is that each microcrack is placed in a non-destructive matrix, but is subjected to an effective stress or strain field, which need not be consistent with the additional far field. Therefore, based on the application of the Mori–Tanaka method proposed by Benvensite [13] for calculating the effective modulus of microcracks, we introduced a microcrack density function as follows:

  Na

2

(17) where N is the number of microcracks per unit area, and a is the half length of the microcracks. Fig. 6 shows a schematic of a model with random microcrack defects. The random microcrack defects were introduced into the DEM model of the SiC ceramics by removing particles to represent the pre-existing random flaw defects. In this model, complex types of defects in the material were simplified to a pre-existing linear flaw. Because the particle ID, which can determine the centre position of the pre-existing flaw, was selected randomly, the distribution of the flaws in the DEM model of the SiC ceramics was random. As shown in Fig. 4, θ is the inclination angle of the flaws measured from the positive direction of the x axis, the value of which could be random ranging from −π/2–π/2. The number and length of the pre-existing flaws are both controllable variables. The number of flaws was assumed N0, and the pre-existing number of flaws N within a unit area was given as follows: N 

N0

(18)

wl

where l and w denote the height and length of the DEM model SiC ceramic, respectively. The length of the pre-existing flaws was given in terms of the product of the particle number n removed from a tip of the flaw to another tip and the minimum particle diameter d. Therefore, the half-length of the pre-existing flaws a was given by the product of the deleted particle number and the minimum particle radius. It can be calculated as follows: (19) a  n  rlo where rlo was the minimum particle radius. Combining Eqs. (17), (18), and (19), the crack-density function can be expressed as follows.

  Na  2

N0 wl

 (n rlo )

2

(20)

Fig. 7 shows the pre-existing flaws in the DEM model of the SiC ceramics. In the model, the numbers of pre-existing flaws N0 were 10, 20, 30, and 40, respectively, and the numbers of deleted particles n were 5, 10, 15, and 20. Table 3 lists the corresponding relationship between the flaw number N0, deleted particle number n, and the crack density α. Table 3

Corresponding relationship between the flaw number N0, deleted particle number n, and crack density α N0 5 10 15 20 n 10 0.002 0.007 0.016 0.028 20 0.004 0.014 0.032 0.056 30 0.005 0.021 0.047 0.084 40 0.007 0.028 0.063 0.113

4.2 Uniaxial compression test

·6·

Fig. 8 shows the uniaxial compression test of the DEM model of the SiC ceramics with the pre-existing random microcrack flaws. In the simulation process, the velocity in the y direction was applied to the specimen by setting a wall at the top and bottom of the DEM model. The wall was represented by a horizontal line with a much larger stiffness than those of the particles to avoid deformation during the loading process. The two walls were moved to the centre of the specimen with a constant v = 0.2 m/s using servo control. Fig. 9 shows the failure mode of the DEM model with the pre-existing flaws under the condition of different crack densities.

5

Results and discussion

5.1 Analysis of failure mode Because of the difference in the number and length of the pre-existing flaws with a random inclination angle, as shown in Fig. 9, the stress distribution in the specimen changed considerably during the loading under the condition of different crack densities. Moreover, this affected the pre-existing random flaws in terms of the crack initiation and propagation, and considerably changed the failure mode of the specimen. When the crack density was small (the crack density is less than 0.028, corresponding to Figs. 9 (a)–(i)), the interaction between the flaws in the pre-existing cracks system is weak. As the uniaxial compression load was increased, new microcracks initiated in the pre-existing flaw tips, continuously merged with the surrounding new initiation cracks or the initiation cracks at the other pre-existing flaw tips, and finally, coalesced into 55º–60º macrocracks. With the increase in the crack density (the crack density was more than 0.028, corresponding to Figs. 9 (j)–(p)), the numerous pre-existing flaws served as stress concentrators and locally transformed the stress state in the model. Because of the influence of the strong stress concentration at the pre-existing flaw tips, the stress in the region with no pre-existing flaw was far less than that at the pre-existing flaw tip. Moreover, the interaction between the flaws gradually increased. After the loading process, the microcracks continued to generate in the flaw tips. However, because of the interaction of the flaws, the generated microcracks spread along adjacent pre-existing flaws at first, and extend approximately in a direction parallel to the maximum principal stress. The failure mode exhibited by the fractured specimen was largely axial splitting. By further observing the failure modes of the specimen with different crack densities, when the value of the crack density approached 0.028, the failure mode of the specimen started to change significantly. To facilitate the analysis of the simulation results, we temporarily considered a crack density of 0.028 as a boundary to divide the crack density into small intervals A (0 < α < 0.028) and large intervals B (0.028 < α < 0.113). Fig. 10 shows the enlarged images of the failure process of the SiC specimen under compression when the crack densities α are 0.007, 0.016, and 0.063. We can clearly observe the entire process from the initiation and growth of the microcracks to the failure of the specimen during the loading. Stages 1, 2, 3, and 4 respectively represent four different loading stages and the corresponding cycles: the loading preparation, initial microcrack, explosive microcrack initiation, and the specimen failure. In the DEM simulation of the uniaxial compression, we assumed that the stress of the specimen decreases below 80% of the peak stress (defined as the compressive strength of the material) as a reference for the end of loading and failure of specimen. Figs. 10a, 10b, and 10c show that the cycles calculated in stages 2 and 4 decreases with the increase in the crack density. The time and velocity of the microcrack initiation become shorter and faster because of the increase in the crack density, and the failure of the specimen become faster and easier. During the loading process, the stress concentration at the crack tips intensified. The microcracks initiated in the pre-existing flaw tip at first. The closer the inclination angle of the pre-existing flaw is to 45º, the easier are the initiation and growth of the cracks (as indicted via the yellow circle in Figs. 10a, 10b, and 10c and the enlarged images of three yellow circles in Fig. 11). In the small

crack-density interval, as shown in Figs. 10a and 10b, the interactions between the pre-existing random flaws was weak. The patterns of the crack initiation, propagation, and coalescence were consistent with those of the centre oblique cracks [27-29]. However, in the large crack-density interval, the cracks generated in the pre-existing flaw tips were affected by adjacent pre-existing flaws. Thus, the expansion patterns of the cracks were largely the propagation and penetration between adjacent pre-existing flaws (Fig. 10c).

5.2 Analysis of stress–strain curve The pre-existing flaws affected the stress distribution of the specimen. Figs. 12 and 13 show the stress–strain curve of the DEM model of the SiC ceramic in the small crack-density interval A and the large crack-density interval B, respectively. During the early loading period, the stress–strain curve was linear, but the different crack densities lead to different slopes in the stress–strain curve. In contrast to the increase in the crack density, the slope of the stress–strain curve was reduced, while the peak stress of the specimen and the corresponding strain values decreased. In the small crack-density interval A, after reaching the peak stress of the specimen, the imposed stress declined dramatically. Consequently, the specimen exhibited brittle failure and lost bearing capacity. However, for the interval B (Fig. 13), increasing the crack density could gradually slow the descending velocity of the curve after reaching the peak stress. In addition, the specimen presented a plastic deformation stage with evident metal-failure feature under compression. Moreover, the volume compression of the specimen transferred to expansion, and the velocity of the axial strain and the volumetric strain increased.

5.3 Analysis of crack initiation number-strain curve Figs. 14 and 15 show the crack initiation number-strain curve of the specimen in intervals A and B, respectively. The microcracks resulted in explosive initiation after reaching the peak stress by constantly loading. The increase in the crack density would promote the explosive microcracks to occur earlier and the specimen failure becomes easier. In the interval A (Fig. 14), as the strain increased to the certain value, the crack initiation curve became so steep that the specimen failure occurred suddenly and the failure mode was brittle fracture. Nevertheless, with the increase in the crack density, the increasing trend of the microcrack initiation number became increasingly smooth in the interval B. At that time, the failure mode of the specimen was similar to plastic deformation. 5.4 Analysis of different random distributed flaws To eliminate the randomness of the established model as much as possible, more DEM tests were performed with five different randomness flaw distributions, as shown in Fig. 16. The failure images of the SiC specimen under compression were presented for crack densities 0.007, 0.016, and 0.063. In the small crack-density interval, as shown in Figs. 16(a) and 16(b), the failure mode of the specimens with five different flaw distributions was brittle fracture principally. However, in the large crack-density interval, as shown in Fig. 16(c), the axial splitting was the main failure mode.

Figs. 17(a), 17(b), and 17(c) show the stress–strain curve of the specimens with five different flaw distributions (distribution 0 is the flaw distribution in Fig. 9) for three different crack densities. For the same crack density, the compressive strength exhibited a slight fluctuation. However, in the three different crack densities, the analysis of the stress–strain curve is consistent with that given in section 5.2. To further accurately study the effect of specimen with random flaws on the macro mechanical properties for different crack densities, we considered the average of the DEM test results of five different flaw distributions in the subsequent analysis (Figs. 18, 19, 20, and 21). The error bar represents the

·8·

standard deviation of the relevant results.

5.5 Analysis of compressive strength The compressive strength directly indicated the mechanical properties of the brittle material in the case of uniaxial compression. Figs. 18 and 19 show the change in the compressive-strength curve and crack-initiation curve as the crack density increased. In Fig. 18, the black and red dots represent respectively the average simulation results of the compressive strength in intervals A and B. The constant increase in the crack density would significantly weaken the compressive strength. Furthermore, its downward tendency in the interval A is clearer than that in the interval B. When the crack density α is changed from 0.002 to 0.028, the compressive strength decreases sharply to only approximately 60% of the original value. As shown in Fig. 19, the downward trend of the crack initiation stress of the specimen in the interval A is even greater than that in the interval B. As the crack density was increased to approximately 0.028, the average crack-initiation stress was reduced by half. The bearing capacity of the specimen reduced significantly by increasing the crack density.

5.6 Analysis of Poisson’s ratio The Poisson’s ratio is defined as the ratio of the transverse normal strain to the axial normal strain under the condition that the material is in the one-way tension or compression. The variation in the Poisson’s ratio with respect to the different crack densities in the DEM simulation was opposite to that of the compressive strength and crack-initiation stress. Fig. 20 shows the Poisson’s ratio for the different crack densities. With the increase in the crack density, the Poisson’s ratio and the transverse strain rate of the material increased. The upward tendency of the Poisson’s ratio of the specimen in the interval A was slightly higher than that in the interval B. In addition, the plastic deformation occurred in the material during the compression process. 5.7 Comparison of effective Young’s modulus The Young’s modulus is the physical quantity that describes the ability of the material to resist deformation, and is defined as the ratio of the stress to the strain of the material. By comparing the effective Young’s moduli with respect to different crack densities obtained from the DEM simulation and the Mori–Tanaka method, Fig. 21 shows the correlation results. The equation to calculate the effective elastic modulus in the Mori–Tanaka method is as follows [14]. E

E



1 1  

(21)

Fig. 21 shows that when the crack density is small (in the interval A), the DEM-simulation curve is strongly consistent with that of the Mori–Tanaka method. However, in the large crack-density interval B, with the increase in the crack density, the effective Young’s modulus simulated using the DEM was slightly lower than the theoretical value of the Mori–Tanaka method. Moreover, the location of the corresponding simulation-fitting curve was below that of the Mori–Tanaka method curve. This is because only the weak interaction of the cracks was considered in the Mori–Tanaka method; thus, this method can ensure sufficient precision only when the distribution of the microcracks is relatively sparse. It can be inferred that the DEM is more precise and is more advantageous for studying materials with defects.

6

Conclusions In this study, the crack-density function wherein the number and length of the pre-existing flaws are considered

variables was introduced to develop the DEM model of the SiC ceramic with pre-existing random flaws based on Mori–

Tanaka method. The effects of random flaws on the failure mode and mechanical properties of the SiC ceramic material were studied. The conclusions of this study are as follows. (1) When the crack density approached 0.028, the failure mode of the DEM model of the SiC ceramic started to change from a brittle fracture to a plastic failure. In the small crack-density interval (the crack density is less than 0.028), the specimen ruptures along angles in the range of 55°–60°. The internal crack extension was relatively independent and the effect of the pre-existing flaws on the same was negligible. The patterns of the crack initiation, propagation, and coalescence were consistent with the centre oblique cracks. In the large crack-density interval (the crack density is more than 0.028); the failure mode was largely axial splitting. The interaction between the cracks was considerable. The generated microcracks spread along adjacent pre-existing flaws at first, and extended approximately in a direction parallel to the maximum principal stress. (2) The pre-existing random flaws have largely weakened the mechanical properties of the materials. With the increase in the crack density, the stress–strain curve of the specimen with a distinct brittle-material characteristic change to that with a metal-failure characteristic under compression. In comparison to the large crack-density interval B, the declining trends of the compressive strength and crack-initiation stress of the specimen in the small-crack density interval A are larger. When the crack density is 0.028, the compressive strength and crack-initiation stress sharply decrease to approximately 60% and 50% of the original value, respectively. (3) For small crack densities, the effective Young’s modulus of the SiC ceramic specimen obtained using the DEM simulation is in good agreement with that obtained using the Mori–Tanaka method. However, with the increase in the crack density, the effective Young’s modulus using the DEM simulation is less than that using the Mori–Tanaka method, wherein only the weak interaction between the cracks is considered.

Acknowledgement This work was supported by the National Natural Science Foundation of China (Grant numbers 51605409, 51605410 and 11602212) and the Provincial Natural Science Foundation of Hunan (14JJ6013).

Reference [1] D. Munz, T. Fett, Ceramics: mechanical properties, failure behaviour, materials selection, Springer Science & Business Media, Berlin, 2013. [2] M. Ortiz, Microcrack coalescence and macroscopic crack growth initiation in brittle solids, Inter. J. Solids. Struct. 24 (1988) 231-250. [3] A. Pandey, A. Shyam, T.R. Watkins, E. Lara-Curzio, R. J. Stafford, K.J. Hemker, The uniaxial tensile response of porous and microcracked ceramic materials, J. Am. Ceram. Soc. 97 (2014) 899-906. [4] M.F. Ashby, S.D. Hallam, The failure of brittle solids containing small cracks under compressive stress states, Acta. Metall. 34 (1986) 497-510. [5] P. Cao, T. Liu, C. Pu, H. Lin, Crack propagation and coalescence of brittle rock-like specimens with pre-existing flaws in compression, Eng. Geol. 187 (2015) 113-121. [6] Kushch. V. I, Microstresses and effective elastic moduli of a solid reinforced by periodically distributed spheroidal particles, Inter. J. Solids. Struct. 34 (1997) 1353-1366. [7] V.I. Kushch, I. Sevostianov, L. Mishnaevsky, Effect of crack orientation statistics on effective stiffness of mircocracked solid, Inter. J. Solids. Struct. 46 (2009) 1574-1588. [8] D. Krajcinovic, Damage mechanics, Elsevier, Holland, 1996. [9] M. Kachanov, Effective elastic properties of cracked solids: critical review of some basic concepts, Appl. Mech. Rev. 45 (1992) 304-335.

·10· [10] Z.P. Bazant, Mechanics of distributed cracking, Appl. Mech. Rev. 39 (1986) 675-705. [11] W. Yang, W.B. Lee, Mesoplasticity and its Applications, Springer Science & Business Media, Berlin, 2013. [12] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta. Metall. 21 (1973) 571-574. [13] Y. Benveniste, On the Mori-Tanaka's method in cracked bodies, Mech. Res. Commun. 13 (1986) 193-201. [14] P.A. Cundall, O.D. Strack, A discrete numerical model for granular assemblies, Geotechnique. 29 (1979) 47-65. [15] H. Huang, B. Lecampion, E. Detournay, Discrete element modeling of tool-rock interaction I: rock cutting, Int. J. Numer. Anal. Met. 37 (2013) 1913-1929. [16] L.A.M. Camones, E.A. Vargas. R.P. Figueiredo, R.Q. Velloso, 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. [17] S. Hentz, F.V. Donze, L. Daudeville, Discrete element modeling of concrete submitted to dynamic loading at high strain rates, Comput. Struct 82 (2004) 2509-2524. [18] M.Y. Zang, Z. Lei, S.F. Wang, Investigation of impact fracture behavior of automobile laminated glass by 3D discrete e lement method, Comput. Mech. 41 (2007) 73-83. [19] W. Gao, M.Y. Zang, The simulation of laminated glass beam impact problem by developing fracture model of spherical DEM, Eng. Anal. Bound. Elem. 42 (2014) 2-7. [20] Y.Q. Tan, D.M. Yang, Y. Sheng, Study of polycrystalline Al 2O3 machining cracks using discrete element method, Int. J. Mach. Tool. Manu, 48 (2008) 975-982. [21] Y.Q. Tan, D.M. Yang, Y. Sheng, Discrete element method (DEM) modeling of fracture and damage in the machining process of polycrystalline SiC, J. Eur. Ceram. Soc. 29 (2009) 1029-1037. [22] Y.Q. Tan, S.Q. Jiang, D.M. Yang, Scratching of Al2O3 under pre-stressing, J. Mater. Process. Tech. 211 (2011) 1217-1223. [23] Y.Q. Tan, S.Q. Jiang, S.J. Nie, Prestress scratching on SiC ceramic, Int. J. Appl. Ceram. Tec. 9 (2012), 322-328. [24] P.A. Cundall, PFC2D user's manual (Version 3.1), Minnesota: Itasca Consulting Group Inc, Minneapolis, 2004. [25] D.M. Yang, Y. Sheng, J.Q. Ye, Discrete element modeling of the microbond test of fiber reinforced composite, Comp. Mater. Sci. 49 (2010) 253-259. [26] D.O. Potyondy, P.A. Cundall, A bonded-particle model for rock, Int. J. Rock. Mech. Min. 41 (2004) 1329-1364. [27] P. Baud, T. Reuschlé, P. Charlez, An improved wing crack model for the deformation and failure of rock in compression, Int. J. Rock. Mech. Min. Pergamon. 33 (1996) 539-542. [28] 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. [29] C.A. Tanga, S.Q. Koub. Crack propagation and coalescence in brittle materials under compression, Eng. Fract. Mech. 61 (1998) 311-324.

List of figure captions: Fig. 1. Microcracks in DEM Fig. 2. Contact model of ball-ball Fig. 3. Parallel bond in PFC Fig. 4. (a) Normal behavior and (b) shear behavior of the bond at contact Fig. 5. Calibration tests of mechanical properties, (a) Uniaxial compression test, (b) single-edge notched beam test and (c) three-point bending test. Fig. 6. Model with random microcrack defects Fig. 7. Pre-existing flaws in DEM model of SiC ceramics Fig. 8. DEM model under uniaxial compression test Fig. 9. Failure mode of specimen under the condition of different crack density (n,N0, α) Fig. 10. Enlarged images of the failure process of the specimen for three different crack densities, (a) 0.007, (b) 0.016

and (c) 0.063. Fig. 11. Enlarged images of (a) yellow circle 'A', (b) yellow circle 'B' and (c) yellow circle 'C' indicated in Fig. 10. Fig. 12. Stress–strain curve in the small crack–density interval Fig. 13. Stress–strain curve in the large crack–density interval Fig. 14. Crack initiation number vs. strain curve in the small crack–density interval Fig. 15. Crack initiation number vs. strain curve in the large crack–density interval Fig. 16. Failure images of specimens for three different crack densities, (a) 0.007, (b) 0.016 and (c) 0.063. Fig. 17. Stress–strain curve of the specimens in five different flaw distributions for three different crack densities, (a) 0.007, (b) 0.016 and (c) 0.063. Fig. 18. Compressive strength vs. crack-density curve Fig. 19. Crack initiation stress vs. crack–density curve Fig. 20. Poission’s ratio vs. crack-density curve Fig. 21. Comparison between DEM simulation and Mori-Tanaka method

·12·

·14·

·16·

·18·

·20·

·22·

·24·

·26·