A numerical investigation of mesh sensitivity for a new three-dimensional fracture model within the combined finite-discrete element method

http://dx.doi.org/10.1016/j.engfracmech.2015.11.006 0013-7944/ 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding author at: Department of Mechanical Engineering, University College London, Torrington Place, London WC1E 7JE, United K Tel.: +44 (0)20 3108 9514. E-mail addresses: liwei.guo@ucl.ac.uk (L. Guo), j.xiang@imperial.ac.uk (J. Xiang), j.p.latham@imperial.ac.uk (J.-P. Latham), b.izzuddin@impe (B. Izzuddin). Liwei Guo a,⇑, Jiansheng Xiang , John-Paul Latham, Bassam Izzuddin b


Introduction
In the field of numerical modelling of fractures in quasi-brittle materials, linear and non-linear elastic fracture mechanics based methods [1][2][3], the extended finite element method (XFEM) [4][5][6] and meshless methods, such as the element free Galerkin method (EFGM) [7,8] have traditionally been in the dominant positions. Due to the discrete nature of fracture and fragmentation behaviour, discontinuum-based numerical methods that are originally used for granular materials, such as the smoothed particle hydrodynamics (SPH) method [9][10][11] and the discrete element method (DEM) [12][13][14] have also become increasingly popular. In actual numerical simulations of engineering applications, the choice of modelling approach should be based on the likely failure mechanism of the material, i.e. whether it is a failure of material, discontinuity or a combination of both [15].
To fully explore and extend the potential of different numerical methods, there is an increasing interest in combining FEM-based and DEM-based methods to converge to a formulation that has the advantage of using the DEM to capture the discrete behaviour during fracture and fragmentation processes while retaining the accurate characterisation of deformation and stress fields using the FEM. It should be noted that the literature mentioned in this section is not meant to be a comprehensive review of numerical methods in fracture modelling, but a tailored one with the focus on using combined FEM and DEM formulations. In this category, different research groups have come up with various strategies in the development http  The three-dimensional fracture model investigated in this paper [22] is a new development in the context of the combined finite-discrete element method (FEMDEM) [23,32]. A simple example of modelling an impact between a fragile sphere (breakable) and a rigid base (unbreakable) is shown in Fig. 1. It can be seen that the sphere is only slightly damaged and still reserves the kinetic energy to bounce off the base when impact velocity is low (10 m/s, Fig. 1a), but breaks into many fragments of different sizes when impact velocity is high (100 m/s, Fig. 1b). In FEMDEM fracture modelling, the entire domain is treated as a multi-body system (e.g. the sphere and the base in Fig. 1) and each discrete body is further discretised by finite element meshes. The FEM formulation is used to simulate continuum behaviour, i.e. the calculation of strains and stresses in finite element domains, and the DEM formulation is used to simulate discontinuum behaviour, i.e. the calculation of contact forces across discontinuities. Comprehensive descriptions of the FEMDEM method can be found elsewhere [23,33]. Regarding three-dimensional fracture modelling using the FEMDEM method, there are three main benefits. Firstly, the interaction between discrete fracture walls can be modelled more realistically and accurately by contact mechanics in DEM algorithms; moreover, other media, e.g. fluid, can be directly introduced between fracture walls for fluid-structure interaction simulations [25]. Secondly, the FEMDEM fracture model can initiate new fractures and furthermore, fracture mechanics energy concepts are used to limit fracture propagation; it has the advantage of not requiring the specification of initial flaws or any predefined fracture patterns (e.g. no initial flaws in the sphere in Fig. 1), which are normally prerequisites for fracture growth models based on linear and non-linear elastic fracture mechanics. Thirdly, due to the addition of the contact detection and interaction algorithms in the DEM formulation, this fracture model is particularly useful when a large number of fragments are generated after impact (e.g. Fig. 1b), and for modelling fracture and fragmentation in multi-body systems [24].
In terms of specific formulations and algorithms, the new three-dimensional fracture model has the following three main features.

1.
A new space discretisation scheme featuring three-dimensional interface elements has been developed. Using this scheme, any arbitrarily shaped three-dimensional domain can be discretised by 4-node tetrahedral elements and 6node interface elements, which are inserted between tetrahedral elements. The material failure criteria are applied to the interface elements whose failure would physically separate tetrahedral elements and generate discrete fracture surfaces. 2. A Mohr-Coulomb failure criterion with a tension cut-off is used to determine the failure state of interface elements. The shear strength is defined as a function of the normal stress acting perpendicular to the shear direction. Therefore, fracturing behaviour in complicated stress fields can be realistically captured.
(b) Impact velocity 100 m/s. 3. The finite element formulation and the discrete element formulation are separated both in the space domain and in the time domain. In the space domain, the FEM formulation is used for continua, while the DEM formulation is used for the interaction across discontinuities, i.e. boundaries of discrete bodies and fracture surfaces. In the time domain, the continuity between tetrahedral elements is only constrained by interface elements before fracture initiation, then after fracture formation the interaction between tetrahedral elements on both sides of the fractures is purely simulated by contact detection and interaction algorithms. The complete governing equation is solved by an explicit time integration scheme.
The detailed algorithms of this three-dimensional fracture model are quite involved and can be found in Guo [22] and Guo et al. [24]. Here only some key concepts related to mesh sensitivity are briefly introduced.
In order to separate tetrahedral elements according to certain failure criteria, a special type of elements -6-node interface elements are inserted between 4-node tetrahedral elements. The deformation in the continuum domain will generate stresses both in 4-node tetrahedral elements and 6-node interface elements. In the interface elements, the stress is calculated from the relative displacement between two triangular faces of adjacent tetrahedral elements. The displacement d and stress r of an interface element are defined as where r is the normal stress component, corresponding to the normal displacement d n , and s is the shear stress component, corresponding to the shear displacement d s .
The three-dimensional fracture model investigated in this paper is similar to the concept proposed by Hillerborg et al. [34], who assumed there exists a plastic zone corresponding to a micro-fractured zone with some remaining ligaments for stress transfer in front of the actual fracture tip. For a single mode I tensile fracture, for example, the transition from the elastic zone to the discrete fracture via the plastic zone is illustrated by Fig. 2. The white area represents the continuum domain that is intact without any fractures, while the discrete fracture is represented by the light yellow area. The orange area is defined as the plastic zone, which corresponds to the displacement range d np $ d nc in the interface elements. At d n ¼ d np , the normal stress r in the interface element reaches its peak value, which is the tensile strength f t in this case.
Ahead of this position (to the right-hand side), the domain is at a strain softening stage (orange area), so the normal stress r decreases from the tensile strength f t to zero at the actual fracture tip, where the normal displacement d n in the interface element reaches its critical value d nc .
Considering the whole transition from the elastic zone to the discrete fracture via the plastic zone, a stress-displacement relation including strain softening effect is used for the interface elements (Fig. 3), which is similar to the combined single and smeared crack model proposed by Munjiza et al. [35]. It should be noted that the normal stress r and shear stress s are calculated following stress-displacement curves of the same shape but different definitions of f, d p and d c on the curves. It is also worth mentioning that the shape of the curve after the peak stress has a very generalised form in Fig. 3, and specific data sets can be used to define the post-peak curve for any quasi-brittle materials. For example, Rougier et al. [21] adopted a different form for their modelling of granite.
The peak stress f in Fig. 3 represents the material strength, so it means the tensile strength f t when calculating the normal stress component r, and shear strength f s for the shear stress component s. In this model, the tensile strength f t is assumed to be a constant, while the shear strength f s is defined by the Mohr-Coulomb criterion with a tension cut-off, where c is the cohesion, / is the internal friction angle, and r n is the normal stress acting perpendicular to the shear direction. Here the engineering mechanics sign convention is used, so tensile stress is positive and compressive stress is negative.
It should be noted that because the normal stress r n cannot exceed the tensile strength f t , the tension cut-off that happens when r n P f t is automatically guaranteed in Eq. (3).
The interface elements will fail when the displacement reaches its critical value d c , which is defined based on the Griffith theory [36]. It assumes that a certain amount of energy is absorbed by the formation of a unit area of the fracture surface in a brittle medium, which is called the fracture energy G f and can be calculated as Therefore, the normal stress r can be calculated following the stress-displacement relation in Fig. 3, where d np is the maximum elastic displacement in the normal direction, d nc is the critical displacement at failure in the normal direction, z is a heuristic softening parameter obtained by curve fitting using experiment data from direct tension tests of concrete [35,37], and can be calculated by Eq. (7).
In actual simulation using this fracture model, the parameters in Eq. (7) are usually chosen as a = 0.63, b = 1.8 and c = 6.0, which are material properties derived from experiment data [38]. D is a parameter defined to quantify the deformation in interface elements, and is given by Eq. (8).
where d sp is the maximum elastic displacement in the shear direction, d sc is the critical displacement at failure in the shear direction. In a similar way, the shear stress s can be calculated by substituting normal displacement d n with shear displacement d s , and other parameters in the normal direction (with subscript n) with the corresponding parameters in the shear direction (with subscript s). For example, different values can be used for the fracture energy G f for tensile and shear modes, but for the numerical tests in this paper (Sections 3 and 4), the shear failure mode is of second order influence on the Brazilian test (indirect tension) and no influence on the pure tension example, so only the fracture energy G f for the tensile mode is used in the simulations. After interface elements fail, discrete fractures will form between tetrahedral elements, using the faces of adjacent tetrahedral elements as fracture surfaces. At this stage, the stress-displacement relation defined in Fig. 3  the failed interface elements. Instead, the contact detection and interaction algorithms in the DEM formulation will be used to simulate the interaction, e.g. normal compression and sliding friction, between fracture surfaces. It should be noted that the mesh size sensitivity of the contact algorithms after fracture formation has been studied elsewhere [39], so in this paper the mesh size sensitivity (Section 3) is investigated by modelling opening-mode fractures and is mainly associated with the FEM formulation. However, in the study of mesh orientation sensitivity (Section 4) the contact algorithms are automatically activated in the modelling of mechanical contact between platens and the disc specimen, and between fracture surfaces in the compressive crushing zones near platens.

Mesh size sensitivity
Previous analytical and numerical studies [31,35] have shown that the size of finite elements close to the fracture tip needs to be much smaller than the length of the plastic zone to achieve accurate results in two-dimensional fracture simulations using the FEMDEM method. In this section, a similar methodology of simulating a series of models with the same geometry but different element sizes is used to investigate mesh size sensitivity for the new three-dimensional fracture model within the FEMDEM method.
In the analysis of numerical results of this section, the length of the plastic zone in front of the actual fracture tip will be used as a main approach to quantify the fracture propagation process. This topic of fracture tip plastic zone in quasi-brittle materials has been extensively studied by analytical and experimental methods [40][41][42][43][44]. In the previous work [31], the theoretical value of the plastic zone length was estimated from analytical solutions [45]. More specifically, the lower value of the plastic zone length D for a short tensile fracture is obtained from Muskhelishvili's solution as where E is the Young's modulus of the continuum, f t is the tensile strength, and d c is the critical opening displacement when bonding stress in the plastic zone equals zero, which can be obtained from Eq. (5), therefore For a long tensile fracture, the lower value of the plastic zone length D is obtained from Westergaard's solution as It should be noted that although Eq. (12) is derived from two-dimensional analysis, the three-dimensional cases simulated in this section can be simplified into two-dimensional problems in theoretical analysis, so Eq. (12) will be used in Section 3.1 as the theoretical estimation for comparison with numerical results.

Test setup
The first problem simulated here is similar to the cases studied in the two-dimensional simulations [31], which is the propagation of a single fracture at the centre of a square domain. The model is shown in Fig. 4. The size of the square domain is 120 mm Â 120 mm in the xy-plane, and 20 mm in the z-direction. A pre-existing horizontal fracture is inserted at the centre of the square. A linearly increasing pressure P is applied at the fracture surfaces and the loading rate is 1:0 Â 10 10 Pa/s. As the pressure increases, the fracture will start to propagate until it breaks the model into two equal parts. As both the geometry and loading condition are symmetric with respect to the central yz-plane, only the right half of the model is simulated (Fig. 4) and a roller boundary condition, which means the translational displacement in the x-direction is constrained, is added to the left boundary of the right-half model.
The material used in the tests (Table 1) is assumed to represent typical rock [46][47][48] or fine concrete mortar properties [49]. The friction coefficient between fracture surfaces is set to be 0.6.
Five models with the same geometry and loading condition but different element sizes are tested. In order to eliminate the influence of mesh orientation, the domain is meshed using structured 4-node tetrahedral elements. Unstructured meshes will be used in Section 4 to investigate the mesh orientation sensitivity. The five meshes are shown in Fig. 5 and the element sizes h and corresponding element numbers N are listed in Table 2.

Numerical results
The numerical results of five models with consecutively refined meshes are compared below. Note that the stress con- Before showing the numerical results, first the theoretical estimations of the plastic zone length are given using the material properties in Table 1. From Eq. (11) the estimation of the lower value of the plastic zone length D lower can be obtained as The upper value of the plastic zone length D upper can be estimated from Eq. (10) Based on the theoretical estimation, for element size of 20 mm (Model 1) and 10 mm (Model 2), the plastic zone can only be discretised by 1-2 finite elements. It can be seen from Figs. 6 and 7 that the numerical results match the theoretical predictions. In these two cases the length of the plastic zone is governed by the element size, which spreads only one element in Model 1 (h = 20 mm) and two elements in Model 2 (h = 10 mm). The stress gradient in front of the fracture tip cannot be accurately captured because there are not enough elements inside the plastic zone.
Refined meshes with element size h = 5 mm (Fig. 8), h = 2.5 mm (Fig. 9) and h = 1.25 mm (Fig. 10) are tested. The results show the length of the plastic zone is independent of the element size. In Model 3 (h = 5 mm, Fig. 8), the plastic zone spreads approximately 3 elements, which is equivalent to 15 mm. The same length of the plastic zone can also be seen in Model 4 (h = 2.5 mm, Fig. 9) and Model 5 (h = 1.25 mm, Fig. 10), which have 6 and 12 elements in the plastic zone, respectively. Especially from the stress contours of element size 2.5 mm (Fig. 9) and 1.25 mm (Fig. 10), the gradient of stress distribution, which decreases from the tensile strength to zero inside the plastic zone, is clearly characterised.
To further compare the plastic zone length obtained from different mesh sizes, the measured values from numerical modelling are plotted in Fig. 11. The normalised element size is the original value divided by 20 mm, i.e. the largest element size  Table 1 Material properties in single fracture propagation tests.

Material properties Values
Density q (kg m À3 )   used in this series of tests. It can be seen that longer plastic zones are generated from larger element sizes and the plastic zone length converges to 15 mm for element size equal to or smaller than 5 mm. Compared with the theoretical estimations given in Eqs. (13) and (14), it can be deduced that the plastic zone should be discretised by at least three finite elements to give a correct numerical representation of the plastic zone ahead of the fracture tip.

Test setup
The second problem is a series of three-point bending tests, i.e. where a beam supported at its two ends is compressed in the middle and in the end it breaks due to flexural deformation. The test setup is shown in Fig. 12. The dimensions of the  beam specimen are 500 mm Â 50 mm Â 20 mm, which are the length, height and thickness in the x; y and z-direction, respectively (Fig. 12a). Loading velocities V = 0.1 m/s are applied in the vertical y-direction to the three platens in order to generate a three-point bending condition. The upper platen moves downwards and the two lower platens move upwards at the same velocity. It should be noted that there is in effect a twofold higher velocity with this setup than a conventional laboratory test where only the central platen moves. To reduce the impact effect, the velocities first increase linearly from zero to a constant value V = 0.1 m/s in 0.1 ms, and then remain constant. The material of the beam specimen in the simulations is assumed to be homogeneous and isotropic, and there are no preexisting flaws inside it (Fig. 12b). The material properties of the beam specimen are the same as the values listed in Table 1. It should be noted that the three-dimensional fracture model is only applied to the beam specimen, and the three steel platens are assumed to be rigid, so material properties are not needed for the platens. The friction coefficient is set to be 0.6 between fracture surfaces, and 0.1 between the beam specimen and platens.  To investigate the influence of mesh sizes on the mechanical behaviour in three-point bending conditions, the beam specimen is discretized by four different mesh sizes (Fig. 13). The element sizes h and corresponding element numbers N are listed in Table 3. All the three platens are meshed in the same manner in four tests, which have 137 elements each.

Numerical results
Fig. 14 shows contours of horizontal stress before fracture initiation and the final fracture pattern in three-point bending tests. From the horizontal stress contours it can be seen that before fracture initiation, the same pattern of stress fields is achieved in all four tests, where the upper part of the beam specimen is in compression and the lower part is in tension, with a neutral surface in the middle of the vertical y-direction. The highest tensile stress happens at the middle of the outer extending arc of the modelled beam, which corresponds to the location of the final fracture (Fig. 14e). It should be noted because all the four tests obtain exactly the same fracture pattern, which is a single fracture breaking the beam specimen into two equal parts, only the final fracture pattern of Beam 4 (h = 2.5 mm) is shown in Fig. 14.      The load F (contact force in the loading y-direction between the upper platen and the beam specimen) is plotted against the maximum deflection d y (vertical displacement at the centre of the beam specimen) for the four mesh sizes in Fig. 15. It should be noted that the contact force is calculated by an integration of nodal contact forces on the platen. It can be seen from the F-d y curves that the peak loads of larger mesh sizes (25 mm and 10 mm) are higher than the peak loads for smaller mesh sizes (5 mm and 2.5 mm), and specimens with coarse meshes fail at smaller deformations. The load-deflection relation and the value of peak load converge to a stable state when mesh size h is equal to or smaller than 5 mm. This observation is in agreement with the results obtained from the same material in single fracture propagation tests (Section 3.1), which showed that the plastic zone length ahead of the fracture tip converges when element size equal to or smaller than 5 mm. It is also in agreement with the conclusion from a two-dimensional numerical simulation of an impact test on a concrete beam with different meshes [35]. It is worth mentioning that despite the over-estimation of peak loads by relatively coarse meshes, the errors are less than 8%, which indicates they might be employed when higher accuracy is not necessary and computational resources are limited.        It can also be seen that the brittle failure behaviour of the beam specimen is correctly captured by the fracture modelling. Here the brittle failure is defined as the significant loss of strength with fracture formation [50]. After reaching the peak value, the load on the beam specimen immediately drops to zero, which means the beam loses its strength to sustain any load so the structure can be regarded as collapsed. It is worth mentioning that there are fluctuations on the F-d y (load-maximum deflection) curves before reaching peak loads; this is because the slim shape of the beam specimen causes certain vibration modes, which affect the recording of the contact force between the beam and the upper platen.

Mesh orientation sensitivity
A complete mesh sensitivity analysis includes two parts: mesh size sensitivity and mesh orientation sensitivity. In the previous section, only mesh size sensitivity is investigated using structured meshes. Once the mesh size satisfies the requirement, the next aspect to consider is mesh orientation. The three-dimensional fracture model used in this paper is based on fixed meshes, so at the element level fractures can only propagate along tetrahedral element boundaries. Tijssens et al. [26] have shown that cohesive zone models show clear mesh dependency of fracture patterns in structured meshes, which means the fractures tend to propagate along dominant directions of element alignment. Therefore, in the three-dimensional fracture modelling using the FEMDEM method, unstructured meshes are recommended in order to reduce the mesh dependency of fracture patterns at the global scale. It should be noted that, even though from a global point of view mesh dependency can be reduced by using unstructured meshes, fracture paths are still dependent on local mesh orientation, and it is necessary to prove the global fracture pattern and critical load are not affected by local mesh orientation when unstructured meshes are used. In this section, specially designed Brazilian tests with unstructured meshes are simulated to examine the mesh orientation sensitivity. It should be noted that the mesh orientation sensitivity studied here does not mean the sensitivity to certain mesh alignment patterns, but refers to the repeatability of numerical results (e.g. fracture path and peak load) using different unstructured meshes with the same mean element size.

Test setup
The setup for the Brazilian tests is shown in Fig. 16. A vertically placed disc specimen perpendicular to the z-direction is compressed diametrically between two platens. The diameter of the disc specimen is 40 mm and the thickness in the zdirection is 15 mm (Fig. 16a). Loading velocities V are applied to both platens to generate an indirect tensile stress field inside the disc. The two loading velocities have the same value but opposite directions. To reduce the impact effect when the loading starts, the velocities first increase linearly from zero to a constant value V = 0.05 m/s in 0.2 ms, and then remain constant in the simulations. The time-step used in the simulations is Dt ¼ 2 Â 10 À9 s.
The domain is meshed using unstructured 4-node tetrahedral elements and the mean mesh size is $1.2 mm. According to the conclusion drawn from Section 3, this mesh size is small enough to generate accurate results so the mesh size effect is not considered in the tests. A total number of 51,690 elements are generated for the disc specimen and 2854 elements for the platens. The two loading platens are originally placed horizontally so the compressive loading is in the vertical y-direction. Then they are rotated with respect to the z-axis to a certain angle h but the disc specimen is kept in its original position (Fig. 16b). The angle between the loading axis and the vertical y-direction is defined as the loading angle h. Because the elements along the loading axis (blue dashed lines in Fig. 16b) are arranged in different patterns when the loading direction changes, the effect of local mesh orientation can be investigated by comparing the fracture patterns and peak loads of four loading angles 0°, 30°, 60°and 90°under identical loading conditions, without the need to actually construct different meshes.
The material properties used for the disc specimen are the same as in the mesh size sensitivity test (see Table 1) except the fracture energy G f is increased to 50 J m À2 . The fracture energy G f was intentionally given a low value in Section 3 because larger fracture energy results in a longer plastic zone (Eqs. (10) and (11)), which is difficult to measure due to the limited size of the domain. It should be noted that the steel platens in the Brazilian tests are assumed to be rigid, so material properties are not needed for them. The friction coefficient is set to be 0.6 between fracture surfaces, and 0.1 between the disc specimen and platens.

Numerical results
The numerical results of loading angles 0°, 30°, 60°and 90°are presented in Fig. 17. It can be seen that although elements are irregularly arranged along the loading axis for different loading angles, the simulations all obtain correct global fracture patterns that match theoretical predictions [51] and the range of experimental observations for homogeneous isotropic rock [52]. Due to the high contact forces, shear fractures first initiate at the two ends of the disc specimen that are in contact with the loading platens. Then the central fracture propagates through the whole disc and splits it into two halves. Final fracture patterns have both major tensile splitting fractures along the loading axis and minor crushing zones (shear fractures) near the loading platens. The fracture path differs somewhat in character in each case and departs more from the diametral loading plane in the h = 90°case (Fig. 17d). However, the failure modes and global fracture patterns are very similar in all the cases, e.g. there are no branches from the middle of the major tensile splitting fracture, which can break the disc into more than two pieces.
The relations between the load F and the axial strain e obtained from numerical simulations are shown in Fig. 18. The load F is calculated as where F 1 and F 2 are the contact forces between the two platens and the disc specimen, respectively. Because both platens move at the same velocity but in opposite directions, the values of F 1 and F 2 are almost equal. The axial strain e is defined to measure the deformation in the disc along the loading axis and is calculated as: where d is the diameter of the disc specimen. From Fig. 18 it can be seen that all the four simulations have the same response at the initial elastic deformation stage.
When the axial strain e exceeds 1.0%, the four curves start to separate and then reach different values of peak loads. The range of peak loads is 2130.0-2266.4 N with a mean value of 2220.0 N, so the variation coefficient (i.e. standard deviation over mean value) is 2.4%, which is comparable for indirect tensile strength of isotropic rock [53]. This shows that if an   unstructured mesh is used and the mean mesh size is small enough, the numerical results of three-dimensional fracture modelling are acceptable regardless of the local mesh orientation.

Computational efficiency
To further test the computational performance of the three-dimensional fracture model, a computational efficiency analysis is conducted using the data recorded from the single fracture propagation tests (see Section 3.1). More specifically, the CPU time is recorded for each of the simulations reported in Section 3.1 and they are compared with respect to the total element number. The current numerical code of the three-dimensional fracture model within the FEMDEM method is a serial code written in C and C++ programming languages. All of the simulations are run on a workstation with Intel Xeon CPU E5-2680 (2.70 GHz). The CPU time needed for one time-step during the fracture propagation stage is plotted in Fig. 19 for five different mesh sizes. It can be seen that the CPU time per time-step increases linearly with increasing element number, which proves that the numerical code works efficiently for different scales. It should be noted that the CPU time needed for FEMDEM modelling also depends on the contact algorithms in the DEM formulation, which might dominate the overall computational performance if there are a large number of discontinuities (e.g. fractures and discrete bodies) in the domain. The research on computational efficiency of the DEM part in the FEMDEM method can be found elsewhere. For example, Munjiza and Andrews [54] studied the contact detection algorithm and reported that the total detection time is proportional to the total number of discrete bodies.   19. Plot of CPU time per time-step versus total element number. The first three data points are also shown in the enlarged window at the right-hand side.

Discussion and conclusions
The mesh sensitivity of a new three-dimensional fracture model within the combined finite-discrete element method was investigated by specially designed numerical tests. Both mesh size and mesh orientation were considered. The sensitivity to mesh size was examined by modelling a single tensile fracture propagation problem and three-point bending tests using a series of models with the same geometry but different structured mesh sizes. The mesh orientation sensitivity was investigated by diametrically compressing a disc specimen from different angles. A very fine and unstructured mesh was used in this test so only local mesh orientation affected the numerical results when the loading angle changed. Moreover, the computational efficiency of the three-dimensional fracture model was studied using the data of CPU time recorded from the mesh size sensitivity test.
From the numerical investigation of mesh size sensitivity it can be demonstrated that the accuracy of three-dimensional fracture modelling depends on the element size around the fracture tips. If the element size is of the same order of magnitude or larger than the theoretical length of the plastic zone, the stress field around a fracture tip is more like a uniform distribution, so the far-field stress has a more significant effect on the fracture propagation than the local stress field. In contrast, for a fine mesh, which can be defined for our purpose as when the element size is only a certain fraction (e.g. one third) of the length of the plastic zone, the gradient of local stress distribution inside the plastic zone can be correctly captured.
The three-dimensional fracture model investigated in this paper is based on fixed meshes and fractures can only propagate along finite element boundaries so the fracture patterns are mesh-dependent. However, the results of the mesh orientation sensitivity test proves that if the element size in an unstructured mesh is smaller than one third of the plastic zone length, although at the element level the fracture path may deviate from the theoretical path to accommodate the element boundaries, from a global point of view, an acceptable solution of the mechanical response of the whole system can still be obtained. Furthermore, if the mesh size is small enough to represent the microstructures (e.g. mineral grains and grain boundaries) in quasi-brittle materials, the roughness of fracture surfaces, rather than being caused purely by mesh dependency can actually represent the realistic microscopic roughness observed in fractured materials.
In general, it can be suggested that unstructured meshes are preferable in fracture simulations for a homogeneous isotropic quasi-brittle material of certain strength properties using the three-dimensional fracture model within the FEMDEM method. Before running an actual simulation, first the theoretical size of the plastic zone should be estimated by Eqs. (10) and (11). Then based on the specific size of the simulated domain, it is essential to choose at least one third of the theoretical plastic zone length as the mean element size in mesh generation. It should be noted that in numerical discretisation of a continuum domain, stress and strain fields in the vicinity of fracture tips are only approximations. In order for the approximation to represent the stress gradient ahead of a fracture tip as accurately as possible, the strategy adopted in this paper is to use low-order (4-node tetrahedron) elements for the whole domain and limit the mesh size around fracture tips. The other approach would be to use high-order finite elements, and in the future it would be worthwhile to explore the possibility of enriching the element libraries to improve the accuracy of the current program. In addition, although the computational time only increases linearly as the total element number grows, large-scale engineering problems may still require unaffordable computational time based on estimates of computational effort from this linear relation. In this respect, parallelisation of the current code combined with the use of less complicated algorithms in non-fractured subdomains might be the most fruitful avenues to provide solutions to overcome this difficulty. Although having those limitations, the current threedimensional fracture model still has the ability to model fracturing and fragmentation behaviour in a wide range of medium-scale engineering problems, such as multi-body collision, fluid-structure interaction and fractured media modelling, in which the whole fracturing process, i.e. pre-peak hardening deformation, post-peak strain softening, transition from continuum to discontinuum, and the explicit interaction between discrete fracture surfaces can be realistically captured.