Certain Discrete Element Methods in Problems of Fracture Mechanics

In this paper two discrete element methods (DEM) are discussed. The free hexagon element method is considered a powerful discrete element method, which is broadly used in mechanics of granular media. It substitutes the methods for solving continuum problems. The great disadvantage of classical DEM, such as the particle flow code (material properties are characterized by spring stiffness), is that they have to be fed with material properties provided from laboratory tests (Young's modulus, Poisson's ratio, etc.). The problem consists in the fact that the material properties of continuum methods (FEM, BEM) are not mutually consistent with DEM. This is why we utilize the principal idea of DEM, but cover the continuum by hexagonal elastic, or elastic-plastic, elements. In order to complete the study, another one DEM is discussed. The second method starts with the classical particle flow code (PFC - which uses dynamic equilibrium), but applies static equilibrium. The second method is called the static particle flow code (SPFC). The numerical experience and comparison numerical with experimental results from scaled models are discussed in forthcoming paper by both authors.


Introduction
The principal problem of classical numerical methods, such as finite element methods, boundary element methods, etc., consists in "too stiff " models, or too complicated simulations of the real states when no a priori knowledge of crack initiation is available. This is why discrete element methods have been introduced to replace fracture mechanics problems by contact problems, which are in many respects more transparent, and which lead us to the same results. Moreover, it is not necessary to know the crack initiation in advance.
In the early 1970s Cundall, [2], and others, [3], introduced discrete elements starting with dynamic equilibrium. First, brick-like elements were used (professional computer program UDEC), and later circular elements in 2D and spherical elements in 3D (PFC -particle flow code -both computer systems issued by ITASCA) simulated the continuum behavior of structures. The application of such methods took place mainly in geotechnics, where soil is a typical grain material with the above-mentioned shape, [11]. If the material parameters are well chosen, the mechanical behavior of the discrete elements is very close to reality. The problem consists of finding such material parameters. There have been many attempts to find these parameters, but still there is no satisfactory output from these studies. A promising approach may be to cover the domain, defining the physical body by hexagonal elements, very close to disks, which can cover the domain with very small geometrical error.
Replacing the discrete elements by elastic, or elastic-plastic, hexagons with the full contact of adjacent elements along their common boundaries yields a honeycomb-like shape of the elements, see, e.g. [15], covering the structure of the continuous medium. It is necessary to note that beams form the honeycomb boundaries in [15] and there is no material inside such particles. In our case some kind of material fills the interior of the hexagonal elements. The relations inside the hexagonal particles are solved by a special form of the boundary element method, [1]. Free hexagons are used by Onck and van der Giessen in [12], for example, where wide range of references on this topic can be found. In [12], the finite element method, e.g. [16], is used to create the stiffness matrices of the elements, namely six finite elements are substructured to a hexagon.
In applications to geotechnical problems the disturbed state concept (DSC) established by Desai, [4], [5], can describe a wide spectrum of material states inside the elements, starting with elastic, elastic-plastic, [6], and even damage states, [10]. The use of eigenparameters for plastic strain, or relaxation stress, [7], [8], completes the description of the possible and suitable nonlinear constitutive laws, which moreover can be "tuned" from "in situ" measurements, or from results from scale modeling. Geotechnical properties are defined on the boundaries of the elements. A typical formulation of the problem involving the generalized Mohr-Coulomb law combined with exclusion of tensile zones is proposed in [14], where the technique using Lagrangian multipliers leads to a mixed problem (both displacements and stresses -element boundary tractions are iterated). In this paper, the penalty method is applied, and element boundary tractions (formerly Lagrangian multipliers) are substituted by spring stiffnesses (i.e., by penalty functions, or in our case by penalty parameters). The springs enable us to simulate the interfacial constraint, namely the exclusion of the tensile tractions and application of the Mohr-Coulomb law. The Mohr-Coulomb law is used in two basic forms, for brittle or almost brittle materials, and for soft rocks or soil.
Several phenomena, e.g., gas extrusion in a coal seam, swelling, watering, and even prestress, can be modeled by Eshelby forces [9]. This treatment seems to be much more promising then the eigenparameters introduced in each element from the zone of some disturbance occurrences, because only tractions (Eshelby forces) along the boundaries of those zones have to be applied.
The paper starts with the formulation of the free hexagon element method, and then statical particle flow is described. Applications in several fields of practical problems are discussed in a forthcoming paper by the authors [13]. tage of some classical DEM, however, is to feed them with material properties provided from laboratory tests (this is the case of the statical particle flow code, formulated in the next section, as the discs are connected by springs, while laboratories provide completely different material parameters). This is here overcome by considering the material characteristics, which are similar to a continuum. The principal idea of classical DEM is adopted, and the domain defining the structure continuum is in our case covered by the hexagonal element, or other material properties can be introduced, such as elastic-plastic, visco-elastic-plastic, etc. This step avoids the necessity to estimate the material properties of springs, which are essential, e.g., for PFC. The free hexagonal element method fulfills a natural requirement due to the fact that the elastic properties are assigned to the particles, and other geotechnical material parameters (angle of internal friction, shear strength or cohesion) to the contacts of the elements. Since most particles are of the same shape it is possible to apply very powerful iteration procedures, because the stiffness matrix can be stored in the internal memory of a computer.
When dealing with crack problems, two principal methods are used. First, the means of fracture mechanics are applied, or secondly a contact problem can be formulated. The first case is generally not suitable, because the direction or way of propagation of the cracks needs to be known in advance.
This paper uses the second possibility, which avoids the obstacle of a priori unknown way of propagation of the cracks by creating a mesh of free hexagonal elements. They are in mutual contact in the undeformed state, but can be disconnected when the contact conditions are violated.
The computational model is described in the following paragraph, where the relations needed for numerical computation are also introduced. The interface conditions are formulated in paragraph 2.2, where the Lagrangian principle is based on the penalty method. The penalty parameters are spring stiffnesses; the springs connect the adjacent elements. The material characteristics of springs can possess a large value to ensure the contact constraints. On the other hand, if, say, the tensile strength condition is violated, the spring parameters tend to zero, and in this case naturally no energy contribution in the normal direction to the element boundary appears in the energy functional in this case. This process excludes the possibility of a multivalue solution, and uniqueness of the solution of the trial problem is ensured. If we cut out the springs when a certain interfacial condition is violated, the problem turns to singular and has not unique solution. Then the way on how the particles move in some later stages of destruction of the trial structure cannot be described.
The hexagonal particles are studied under various contact (interfacial) conditions of the grain particles (elements). In our paper two contact conditions are considered: l the generalized Mohr-Coulomb hypothesis, with exclusion of non-admissible tensile stresses along the contact (a rock mass), l limit state of shear stresses and exclusion of tensile tractions along the contact (a brittle coal seam).
The first case is generally connected with applications in geotechnics, composite materials, shotcrete, etc., and the second case is more appropriate for applications in underground bumps or rock bursts. A two-dimensional formulation and its solution have been prepared and are studied in this paper.
The problem formulated in terms of hexagonal elements (which are not necessarily mutually connected during the loading process of the body, because of nonlinearities arising due to the interfacial conditions) enables us to simulate that crack propagate. The cracking of the medium can be described in such a way that the local damage may be derived. Local deterioration of the material is also shown in the pictures drawn for particular examples. Such a movement of elements and change of stresses probably cannot be obtained from continuous numerical methods.

Computational model
Let us now consider a single hexagonal element (described by domain with its boundary ). Its connection with the adjacent elements is shown in Fig. 1. In each hexagonal element, the pseudo-elastic material properties are taken into consideration, i.e., in each iteration step the element behaves linearly, but the material properties can change during the process of loading and unloading. This makes it possible to introduce only an elastic material stiffness matrix, which is homogeneous and isotropic, and we get well-known integral equations that are valid along the boundary abscissas of the hexagons, [1]: where b i are components of the volume weight vector, G s are edges (abscissas) of the boundary elements, x is the point of observation, x is the integration point, u i are components of the vector of displacements (defined not exclusively on the boundary, but also in the domain of the hexagonal element), p i are components of the tractions, c kl are components of the matrix, the values of which depend on position of the point of observation. The quantities with an asterisk are the given kernels. The kernels can be expressed as (for example, see [1]): , , Assuming uniform distribution of the boundary quantities (displacements u i (x) and tractions p i (x), i = 1, 2, and volume weight forces b i to be uniform in the domain W, and positioning the points of observation x successively at the points x s , which are the centers of the boundary abscissas of the hexagonal elements, a simplified version of (1) is written as: Moreover, the vector of influences of the volume weight forces on the boundary abscissas is , , s = 1, …, 6, and For better and more convenient computation, the most important integrals are given in the Appendix. In this way, the integrals in (2) may be calculated directly, without numerical integration.
Let us introduce vectors a s , b s , s = 1, …, 6, and also u and p as: , u Using this notation, the relations on the elements (2) can be recorded as: where A and B are (12 * 12) matrices, and their components are singular integrals over the boundary abscissas. Matrix A is generally singular, while matrix B is regular. This fact enables us to rearrange equations (3) into the form: where the stiffness matrix K is different from that arising in applications of finite elements (here it is prevailingly non-symmetric), V is the vector of volume weight forces concentrated on the boundary abscissas (more precisely at the point x). In this way, the discretized problem becomes a problem similar to the FEM. Along the adjacent boundary abscissas it should hold (p i m are Eshelbys' forces): where superscript plus means from the right and minus from the left (at most two particles can be in contact). Now using the relations (4) and (5), we get twice as many unknowns as equations, because no connection between the elements has yet been introduced. Equations (5) have to be accomplished by a constraint of the type The latter conditions are penalty-like conditions, since if k i is great enough, the distribution of displacements is continuous, and the displacement from the right is equal to the displacement from the left. These conditions can locally be violated, because of the contact conditions, which are discussed later in this text. Introducing boundary conditions and assuming that k i remains great enough leads us to a stable system of equations delivering a unique solution. Even in the case when local disturbances occur, the solution can be stable. It can happen that there are too many disturbances, e.g., dense occurrence of crack, and localized damage along a path (earth slope stability violation). Then the solution is unstable, and there is a failure of the structure. This is also, for example, the case of a rock burst.
Discretization in the previous sense leads to a nonlinear system of algebraic equations, which are solved by an over-relaxation iterative procedure. This method is sufficient for study purposes. For a larger range of equations the conjugate gradient method has been prepared.
For displacements inside the element domain W it holds: where the element boundary displacements and tractions are known from the previous computation, providing the solution is stable. Using kinematic equations and Hooke's law, the internal stresses can be calculated from (7). There is no danger of singularities, as the points x and x never meet (point x lies inside the domain W and x on boundary G).

Formulation of the contact problem
Recall that displacements are described by a vector function ( ) u = u u 1 2 , of the variable ( ) x = x x 1 2 , . The traction field on the particle boundaries is denoted either as ( ) p = p p 1 2 , , or after projections to normal and tangential directions as . A similar result is valid for projections of displacements, .
Assuming the "small deformation" theory, the essential contact conditions on the interface may be formulated as follows (no penetration conditions): where G C k , k = 1, …, n are boundaries between adjacent particles, u n k, a is the normal displacement of current element a = c and a = a belongs to the adjacent element, both on the current common boundary , G C k , k runs numbers of all common sides of the particles, n is the number of common sides of hexagons (having exactly two adjacent particles inside the domain, one or none on the external boundary).
Let k n k be the spring stiffness in the normal direction and k n k be the spring stiffness in the tangential direction on the boundary between particles with a common boundary G C k .
Then in the elastic region  is the set of displacements that fulfill the kinematical boundary conditions and condition (8). If p n k = 0 then set K is a cone of admissible displacements satisfying the essential boundary and contact conditions. This is valid for brittle or almost brittle material. If the material exhibits elastoplastic behavior, then cone K is changed as: where f is the angle of internal friction, and p n k is the normal traction on the side k, c is the generalized Heaviside's function being equal to zero for a positive argument and equal to one otherwise. Here the sign convention is important: positive normal traction is tension.
where e is the strain tensor, C is the stiffness matrix of the particle, T denotes transposition, W 0 is the sum of subdomains W, i.e., of hexagonal elements, b is the volume weight vector.
Note that the spring stiffness k n k plays the role of a penalty. Recall that the problem can also be formulated in terms of Lagrangian multipliers, and then leads to mixed formulation. The latter case is more suitable for a small number of boundary variables; the problem discussed in this paper decreases the number of unknowns introducing the penalty parameters.

Statical PFC
This section deals with the idea of modeling a structure, e.g., an earth body, by using a statical version of PFC. Recall that the PFC is based on dynamic equilibrium; for slow movements of structures, which appear in most geotechnical problems, for example, it seems to be better to employ statical equilibrium. The earth mass is modeled in a statical version by balls in 3D or disks in 2D in a similar way as in the dynamic version. The balls are connected by springs that relate forces and the appropriate difference of displacements in the direction of the springs. The springs are considered either in normal and tangential directions to the boundary of the particles, or in the direction of the coordinate axes x and y.
In our case, statical equilibrium has to be fulfilled all over each ball, and at the contact points between the adjacent balls. The balls (disks) are considered rigid. Introduce coordinate system 0xyz in 3D. Then each disk has six degrees of freedom (displacement u in direction x, displacement v in direction y , and displacement w in direction z, and three rotations with respect to the three axes x, y and z. In what follows we restrict our considerations to 2D for simplicity; generalization to 3D is straightforward. The movement of each disk is described by two displacements u, v and rotation j. The forces concentrated at one contact point of the adjacent balls obey the contact conditions that are typical for soil in our study. They determine the change of stages in the DSC model. The plastic behavior, providing rigid balls, is imposed only by the forces brought about by spring stiffnesses and eigenparameters, e.g., plastic strains (displacements), or relaxation stresses (forces). When introducing some spring (of the shape of a straight line, or more precisely of an abscissa) spanned between two points, its stiffness is k, the relation force F -displacement u in the elastic case may generally be written as (considering one-dimensional case): where respectively l is the eigenstress (eigenforce), and m is the eigenstrain (eigendisplacement). Another possibility is to drop out the eigenparameters and impose the nonlinear conditions exceptionally on the springs. The eigenparameters enable a larger range of physical laws to be used in the material models. This is why they are mostly applied in the mechanical models.

Relations between disks
Let us take a set of disks describing a discontinuous medium positioned in a coordinate system 0xy, see Fig. 2. At the contact points (nodal points), the springs in the tangential and normal directions are introduced, which have stiffnesses k t in the tangential direction to the boundary of the disk and k n in the normal direction to the boundary at the interface nodal point connecting two adjacent disks. Such a connection is described in more detail in Fig. 3 for three disks in mutual contact. Fig. 3 depicts external (volume weight) forces F i , i = 1,…, n (n is the number of all disks), contact forces in normal and tangential directions N ij , Q ij , respectively, where i, j = 1, 2, and reactions at the supports (created in this case by a flat plate in 3D, or by a straight line in 2D) A, B, H A , H B . The first index denotes the number of the current disk and the second index stands for the number of the adjacent disk. This notation is kept in the following text. The main objective is to formulate the equations of equilibrium in each disk i, i = 1, …, n and from this equilibrium to determine the displacements of the center u i , v i and the rotations j i of each disk. The connection with the adjacent disks is created by the quantities with indices i (the current disk) and j (the adjacent disk).
In the sense of (10)  A typical disk with springs introduced at nodes in x and in y directions and the induced forces in the normal and tangential directions are illustrated in Fig. 4.
If no rotations were considered, the above formulas would be valid without improvement and the computation may start with (13). For each disk two degrees of freedom in 2D (two independent displacements are unknown), and three DOF in 3D (three independent displacements) are sought.
In the case of admitted rotations of disks, additional unknown angles j i describing the rotations of the disks have to be introduced. Recall that three DOF (two displacements  (15) where w ij is given for each node and r i is the radius of disk i.
The forces N x ij , and N y ij acting between disk i and disk j are then given by (13).

Governing equations
Equation (13) can be rearranged in a more suitable way for algorithmization: The third unknown j i appears in the conditions of equilibrium in nonlinear terms, namely in cosines and sines (recall that w ij is the angle of deviation from the x-axis of the point ij on the boundary of disk i being in contact with disk j). In order to avoid very complicated and unreliable nonlinear computation, the load (for example volume weight) will be divided into increments, and in each increment the small displacement (or more precisely small rotation) theory will be considered. Assuming small enough increments, small enough angle j also results, and (17) is substituted by: the possibility to test the contact conditions, as described in the following section. An additional condition is necessary to complete the equations for three unknowns in each disk. This is the moment condition with respect, for example, to the center of the disk under study: where M i is the number of nodes on the boundary of disk i.
Since r i is constant in each disk i, it may be dropped out and the three conditions in disk i are obtained as: where F x i and F y i are volume weight forces. If only gravitation is considered, the denotation F F i y i º , and assumption F x i º 0 can be used, as in Fig. 3.
After introducing boundary conditions and eigenforces, the equations (19) create an algebraic system for unknown displacements u x,tran i , u y,tran i , and rotations j i of the disks.
Solving this, the forces can be determined from (16). Note that when properly stated, system (19) has a unique elastic solution, providing that the incremental formulation (17a) is assumed, i.e., small displacement theory may be employed.

Interfacial conditions
In geotechnics and mining engineering it is well known that the material behavior of the soil mass is not elastic, but exhibits either plastic behavior, or localized damage, or both. Contrary to the case of free hexagons, plastic behavior has to be employed only for spring properties. For the sake of completeness, an example of such properties is discussed here.
Consider two balls that are connected by the above-described system of springs. For the sake of simplicity, we again concentrate our attention on the 2D case. Using formula (11), the normal forces N ij and tangential forces Q ij can be obtained for each linear state. Vector transformation of the coordinates applied to (17) provides displacements u n ij and u t ij .
In a wide range of problems, the classical plastic laws of deformation and stress state, e.g., elasto-plastic law, are formulated as: where s is an exponent to be stated from laboratory tests, as well as the starting value of k n0 ij . The value of N ij increases nonlinearly and when the strength is reached, a crack is assumed at point ij and the spring is suddenly removed.
In practical examples, the spring that is in tension is removed gradually to stabilize the convergence and speed up the iteration process. On the other hand, at this point "penetration" of one disk into the adjacent disk is fully permitted. This is an impact of the nonlinear behavior (21) of the spring being compressed in the normal direction. Damage occurs when the Mohr-Coulomb hypothesis is violated or tensile strength is reached. In this case, it means that a) ( ) where t b is the shear strength (cohesion), c is the Heaviside function, and f is the angle of internal friction, where N ij + is the tensile strength.
When condition a) is not fulfilled, then a "cut" of Q ij is supposed according the following rule: Note that more complicated rules may be imposed. For example, both internal parameters, angle of internal friction and shear strength, may change with the values of Q ij . In the case of violation of condition b), a local disconnection (debond) occurs and the spring is removed again. This is not due to local cracking but because of disconnection of the disks that were originally in contact at point ij.
Eigenforces have not been discussed in this paper. They are additional design parameters for optimal approximation of reality or of laboratory tests in such a numerical model. They may also simulate other phenomena, such as change of temperature, gas emission, blasting at some local sources, etc.

Conclusions
Two new discrete element methods have been introduced in this paper. One of them is an extension of the well--known particle flow code. The solution of this method is very easy and fast. On the other hand, it bears the same disadvantages as PFC itself. The interpretation of the material properties currently obtained from laboratory tests is quite complicated, if possible at all. If the material properties are properly chosen, the results seem to be realistic, as shown in the forthcoming paper by the authors. A much more complex and suitable method for predicting the real behavior of the test material is the method of free hexagons, which involves both mechanical and geotechnical properties received from the experiments. The hexagons can cover the entire domain describing the physical body, and along the local boundaries between elements the geotechnical properties can be imposed. The material behavior inside the elements is described by virtue of boundary elements, which are more appropriate than finite elements in this case. When using the finite element method, the local tractions are polynomials of one order higher than the boundary displacements, while the boundary element method delivers polynomials of the same order along the boundary. The spring condition is then better fulfilled by the boundary elements.