AP02_04.vp 1 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 simula- tions 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 trans- parent, 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], intro- duced 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 contin- uum 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 mate- rial parameters are well chosen, the mechanical behavior of the discrete elements is very close to reality. The problem con- sists 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-plas- tic, 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 con- tinuous 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, start- ing with elastic, elastic-plastic, [6], and even damage states, [10]. The use of eigenparameters for plastic strain, or relax- ation 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 com- bined 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 para- meters). 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 ele- ment 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 dis- cussed in a forthcoming paper by the authors [13]. 2 Free hexagonal element method The discrete free hexagonal element method may be considered a discrete element method. The great disadvan- 42 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 42 No. 4/2002 Certain Discrete Element Methods in Problems of Fracture Mechanics P. P. Procházka, M. G. Kugblenu 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. Keywords: Discrete element methods, free hexagon element method, statical particle flow code. 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 laborato- ries provide completely different material parameters). This is here overcome by considering the material characteristics, which are similar to a continuum. The principal idea of classi- cal DEM is adopted, and the domain defining the structure continuum is in our case covered by the hexagonal ele- ment, 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 ap- ply 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 discon- nected when the contact conditions are violated. The computational model is described in the following paragraph, where the relations needed for numerical compu- tation 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 ex- cludes 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: � the generalized Mohr-Coulomb hypothesis, with exclusion of non-admissible tensile stresses along the contact (a rock mass), � 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 under- ground 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 de- scribed in such a way that the local damage may be derived. Local deterioration of the material is also shown in the pic- tures drawn for particular examples. Such a movement of elements and change of stresses probably cannot be obtained from continuous numerical methods. 2.1 Computational model Let us now consider a single hexagonal element (de- scribed by domain with its boundary ). Its connection with the adjacent elements is shown in Fig. 1. In each hexagonal ele- ment, 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]: � � � � � � � � � � � � c u p u u p l ss kl l i ik i id d � � � � � � �� � � � � � 1 2 x x x x x x* *, , �� � � � � � � � � � � � �� � �� �� is i b u k 1 2 1 6 1 2 1i ik dx x x * , ,� � , ,2 (1) where bi are components of the volume weight vector, �s are edges (abscissas) of the boundary elements, � is the point of © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 43 Acta Polytechnica Vol. 42 No. 4/2002 Fig. 1: Geometry of adjacent hexagonal elements observation, x is the integration point, ui are components of the vector of displacements (defined not exclusively on the boundary, but also in the domain of the hexagonal element), pi are components of the tractions, ckl 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]): � � u A M r x x p A k n x n x ik * ik i k 2 ik * 2 k i i k r r � �� � � � � � � � � log , 2 k x x x n�ik i k 2 j jr �� � � � � � 2 where � � � � � � � �A M� � � � � � �� � �� � � � � � �4 2 3, , � �k x x r x x� � � � � �� � � �, , ,i i i 2 1 2 2 2 and � and � are Lame’s material constants. Assuming uniform distribution of the boundary quantities (displacements ui(x) and tractions pi(x), i 1, 2, and volume weight forces bi to be uniform in the domain �, and position- ing the points of observation � successively at the points �s, which are the centers of the boundary abscissas of the hexago- nal elements, a simplified version of (1) is written as: � � � � 1 2 1 2 1 u p u u p ss is k t i s ik i s ikd d� � � ��� �� * *, ,x x x x� � �� � � 6 1 2 1 2 1 6 � �� � � � � � � � b u k t i i s ik d * , , , ; , , ,x x� � � (2) where ui s and pi s are the values of the relevant quantities posi- tioned at the �s, s 1, …, 6, i.e., � �u ui s i s� � and � �p pi s i s� � . Moreover, the vector of influences of the volume weight forces on the boundary abscissas is � �bs � 1 2, , s 1, …, 6, and � � � � k s i s ik s d� ��� � b u k i * , , ,x x� � �1 2 1 2. 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 numeri- cal integration. Let us introduce vectors �s, �s, s = 1, …, 6, and also u and p as: � � � � � � � � s s s s s s� � � � � � � � � u u p p 1 2 1 2 1 2 3 4 5 6 , , u � � � � � � � , p � � � � � � 1 2 3 4 5 6 � � � � , b b b b b b b 1 2 3 4 5 6 � � � �� �k s i s ik s k s i s ik s s s d d� � ��� � � � u p p u i i * *, , ,x x x x� � � �1 2 1 2 � . Using this notation, the relations on the elements (2) can be recorded as: Au Bp b� � (3) 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: Ku p V K B A V B b� � � �� �, ,1 1 (4) 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 �). In this way, the discretized problem becomes a problem similar to the FEM. Along the adjacent boundary abscissas it should hold (pi � are Eshelbys’ forces): � � � �p p p pi i i i� � � � � � �� � , (5) 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 � �k u u pi i i i� �� � . (6) The latter conditions are penalty-like conditions, since if ki is great enough, the distribution of displacements is continu- ous, 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 dis- cussed later in this text. Introducing boundary conditions and assuming that ki 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-re- laxation 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 � it holds: � � � �u p u u p is k i s ik * i s ik * d s � �� � � � � ��� �� x x x , �1 2 1 6 � � � �, , , , ,� �d d s i s ik *x x x � � � �� � � � � � � � b u k i 1 2 1 2 (7) where the element boundary displacements and tractions are known from the previous computation, providing the solu- tion is stable. Using kinematic equations and Hooke’s law, the internal stresses can be calculated from (7). There is no 44 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 42 No. 4/2002 danger of singularities, as the points x and � never meet (point x lies inside the domain � and � on boundary �). 2.2 Formulation of the contact problem Recall that displacements are described by a vector func- tion � �u � u u1 2, of the variable � �x � x x1 2, . The traction field on the particle boundaries is denoted either as � �p � p p1 2, , or after projections to normal and tangential directions as � �p � p pn t, . A similar result is valid for projections of dis- placements, � �u � u un t, . Assuming the “small deformation” theory, the essential contact conditions on the interface may be formulated as follows (no penetration conditions): � �u n k n k, c n k, � � �u u � 0 on �C k , (8) where �C k, k 1, …, n are boundaries between adjacent parti- cles, un k, � is the normal displacement of current element � c and � a belongs to the adjacent element, both on the current common boundary , �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 in- side the domain, one or none on the external boundary). Let kn k be the spring stiffness in the normal direction and kn k be the spring stiffness in the tangential direction on the boundary between particles with a common boundary �C k. Then in the elastic region � �p kn k n k n k� u and � �p kt k t k t k� u . Denote � � � � � � � � K V p p k p p p k � � � �� � � u u u , ,n + k n k n k n k n + k n k n k t k t k if then 0 � � � �� � � � � � �c k n u uk C k n + k t k t k , c t k ,on � , , , , .1 u � where ut k, � is the tangential displacement on the side k, � �pn+ kdenotes the tensile strength, ck is the shear strength, V is the set of displacements that fulfill the kinematical bound- ary conditions and condition (8). If pn 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: � � � � � � � � K V p p k p p p k � � � �� � � u u u , ,n + k n k n k n k n + k n k n k t k t k if then 0 � � � � � � � � � � � � � � c p p k n u u k n k n k C k n + k t k t k , c t k , � � tan on � , , , ,1 u . where � is the angle of internal friction, and pn k is the normal traction on the side k, is the generalized Heaviside’s func- tion being equal to zero for a positive argument and equal to one otherwise. Here the sign convention is important: posi- tive normal traction is tension. From the above defined spaces we can deduce that � �pn k n k, u , and � �pt k t k, u behave linearly between certain limits, which are given by the material nature of the body. The total energy J of the system reads: � � � � � �� � � � J a k a C k n u u u u b u u u � � � � � � � �� � � 1 2 2 1 , , n k n k d d C T T � � � � � � d� � � 0 1 2� � � � � � � , , ,� � � � � � � � � u x u x u x u x n n t t t n n t � � . (9) where � is the strain tensor, C is the stiffness matrix of the particle, T denotes transposition, �0 is the sum of sub- domains �, i.e., of hexagonal elements, b is the volume weight vector. Note that the spring stiffness kn 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 for- mulation. 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. 3 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 move- ments of structures, which appear in most geotechnical prob- lems, 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 direc- tion 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 �. The forces concentrated at one contact point of the adja- cent 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 im- posed 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): F k u� � �, or � �F k u� � � (10) where respectively � is the eigenstress (eigenforce), and � is the eigenstrain (eigendisplacement). Another possibility is to drop out the eigenparameters and impose the nonlinear con- ditions exceptionally on the springs. The eigenparameters © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 45 Acta Polytechnica Vol. 42 No. 4/2002 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. 3.1 Relations between disks Let us take a set of disks describing a discontinuous me- dium 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 kt in the tangential direction to the boundary of the disk and kn 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 Fi, i 1,…, n (n is the number of all disks), contact forces in normal and tangential directions Nij, Qij, 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, HA, HB. 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 equilib- rium in each disk i, i 1, …, n and from this equilibrium to de- termine the displacements of the center ui, vi and the rotations �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) the physical equations at each nodal point are formulated as: N Q k k ij ij n ij t ij n ij t ij n0 0 � � � � � � � � � � � � � � � � � ij t ij � � � � (11) where �n ij and �t ij are respectively eigenforces in normal and tangential directions. Indices i, j describe the numbers of disks in mutual contact and �n ij , �t ij are the differences of dis- placements in normal and tangential directions, respectively, between disk i and disk j, i.e., �n ij n ij n ji� �u u , �t ij t ij t ji� �u u . Let the nodal point ij under consideration be deviated from the x-direction by an angle �ij. Then the transformation of forces to the 0xy coordinate system is written by: N N Nx ij y ij ij ij ij ij i� � � � �� � � � � � cos sin sin cos � � � � j ij ij T ij ijQ N Q � � � � � � � T (12) where Nx ij, Ny ij are forces in x and y directions, Tij is the trans- formation matrix and superscript T denotes transposition. Recall that matrix Tij is unitary, which means that T Tij ij T� �1 . Since the same equations hold for displacements, the following forces – displacements relation holds valid as: N N k k x ij y ij ij T n ij t ij n ij t ij 0 0 � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � n ij t ij xx ij xy ij yx ij yy ij k k k k � � n ij t ij n ij t ij � � � � � � � � � (13) where k k kxx ij n ij ij t ij ij� �cos sin2 2� � ; k k kyy ij t ij ij n ij ij� �cos sin2 2� � ; � �k k k kxyij yxij nij tij ij� � �12 2sin � � � � � �x ij n ij ij t ij ij� �cos sin ; � � � � �y ij n ij ij t ij ij� �sin cos (14) and �x ij , �y ij are differences of displacements in x and y directions, respectively, in disks i and j, i.e., �x ij x ij x ji� �u u , �y ij y ij y ji� �u u . 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 �i describing the rotations of the disks have to be introduced. Recall that three DOF (two displacements ui, vi and one angle of rotation �i in 2D are to be sought. The situation is depicted in Fig. 5. Let us focus on one typical disk. Its basic movements and their denotations are clear from Fig. 6. With respect to the 46 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 42 No. 4/2002 Fig. 2: Structure of the balls Fig. 3: Forces between disks above-mentioned arguments, the three movements and their parts are described in this picture. The total displacement at any point on the boundary of the disk in the x direction, or in the y direction, consists of two kinds of displacements that are related to the rota- tion (subscript rot), and translation (subscript tran), i.e., u u ux i x, rot i x, tran i� � , u u uy i y, rot i y, tran i� � . Since the 2D case is under consideration, the unknown quantities in each disk i are ux, tran i , uy, tran i and �i. It only remains to express the influence of rotation of each node. From Fig. 6 it obviously holds: � �� � � � u r u r x, rot i y, rot i � � � � � � � i ij i ij i ij i cos cos , sin � � � � �� �sin ,�ij (15) where �ij is given for each node and ri is the radius of disk i. The forces Nx ij, and Ny ij acting between disk i and disk j are then given by (13). 3.2 Governing equations Equation (13) can be rearranged in a more suitable way for algorithmization: N N N N k k k k x ij y ij x ji y ji xx ij xy ij xx ij x� � � � � � y ij xy ij yy ij xy ij yy ij xx ji xy ji xx ji xy ji x k k k k k k k k k � � � � � y ji yy ji xy ji yy ji x i y i x j y j � � � � � � � � � � � � � � � � k k k u u u u � � � � � � � � � � x i y i x j y j (16) and from (15): u u u u u u u x i y i x j y j x,tran i y,tran i x,tran � � � � � �� � j y,tran j i ij i ij i u r r � � � � � � �cos cos sin � � � � �� � � �� � � �� � � � � � � � � � � ij i ij j ji i ji j ji i ji � � � � � � � sin cos cos sin sin r r � � � (17) The third unknown �i appears in the conditions of equi- librium in nonlinear terms, namely in cosines and sines (recall that �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 � also results, and (17) is substituted by: u u u u u u u x i y i x j y j x,tran i y,tran i x,tran � � � � j y,tran j i i ij i i ij j i u r r r � � � � � � � � � sin cos sin cos � � � ji j i jir � � � (17a) In each increment ux,tran i , uy,tran i , and �i are parts of the total values, as also are the forces computed from these movements. Another advantage of this incremental process is © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 47 Acta Polytechnica Vol. 42 No. 4/2002 Fig. 4: Starting model of one disk for further computation Fig. 5: Additional displacements at the point ij from rotation Fig. 6: Three basic movements of a disk the possibility to test the contact conditions, as described in the following section. An additional condition is necessary to complete the equa- tions for three unknowns in each disk. This is the moment condition with respect, for example, to the center of the disk under study: Q ij i i r j M � � � 1 0 , (18) where Mi is the number of nodes on the boundary of disk i. Since ri is constant in each disk i, it may be dropped out and the three conditions in disk i are obtained as: N F N F Q N j M j M j M x ij x i y ij y i t ij x ij i i i � � � � � � � � � � 1 1 1 , , sin �� �ij yij ij i � � � � N j M cos � 1 0 (19) where Fx i and Fy i are volume weight forces. If only gravitat- ion is considered, the denotation F Fi y i� , and assumption Fx 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 ux,tran i , uy,tran i , and rotations �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. 3.3 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-de- scribed system of springs. For the sake of simplicity, we again concentrate our attention on the 2D case. Using formula (11), the normal forces Nij and tangential forces Qij can be obtained for each linear state. Vector transformation of the coordinates applied to (17) provides displacements un ij and ut ij. In a wide range of problems, the classical plastic laws of deformation and stress state, e.g., elasto-plastic law, are formulated as: N Q kij ij 0 2 2 24� � , i.e., � � � �k u k u knij nij tij tij 0 2 2 24� � (20) where k0 is a given positive number. In the case of violating conditions (20), new kt ij may be determined. While the value of kt ij is restricted by (20), un ij is mostly bounded by some value that cannot be exceeded. Moreover, kn ij changes nonlinearly with un ij. A parabolic rule is mostly applied as: � �k k unij n0ij nij s � �1 (21) where s is an exponent to be stated from laboratory tests, as well as the starting value of kn0 ij . The value of Nij 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 re- moved 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 per- mitted. 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) � �Q N Nij ij b ij� � �tan � � , where �b is the shear strength (cohesion), is the Heaviside function, and � is the angle of internal friction, b) N Nij ij +� , where Nij + is the tensile strength. When condition a) is not fulfilled, then a “cut” of Qij is supposed according the following rule: � �Q N Nij ij b ij tij� � � �tan � � � . Note that more complicated rules may be imposed. For exam- ple, both internal parameters, angle of internal friction and shear strength, may change with the values of Qij. In the case of violation of condition b), a local discon- nection (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 dis- advantages 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 im- posed. 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 ele- ment method, the local tractions are polynomials of one 48 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 42 No. 4/2002 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. Acknowledgment This research was supported by the Grant Agency of the Czech Republic – grant number 103/00/0530 and by the Ministry of Education of the Czech Republic, project MSM:210000001,3. References [1] Brebbia, C. A., Telles, J. F. C., Wrobel, L. C.: Boundary Element Techniques. Berlin: Springer-Verlag, 1984. [2] Cundall, P. A.: A computer model for simulation of progressive large scale movements of blocky rock systems. Symposium of the international society of rock mechanics, 1971, p. 132–150. [3] Cundall, P. A., Strack, O. D. L.: A discrete numerical model for granular assemblies. Geotechnique, 1979, p. 47–65. [4] Desai, C. S.: A Consistent Finite Element Technique for Work-Softening Behavior. In: J. T. Oden et al (eds.), Int. Conf. on Comp. Meth. in Nonlinear Mech. Univ. of Texas at Austin, 1974. [5] Desai, C. S.: Constitutive Modeling Using the Disturbed State Concept. Chapter 8, Continuum Models for Materials with Microstructure, ed. H. Mulhaus, UK: John Wiley & Sons, 1994. [6] Duvaut, G., Lions, J.-L.: Les inequations en mecanique et en physique. Paris: Dunod, 1972. [7] Dvorak, G. J., Procházka, P. P.: Thick-walled Composite Cylinders with Optimal Fiber Prestress. Composites, Part B, 27B, 1996, p. 643–649. [8] Dvorak, G. J., Procházka, P. P., Srinivas, S.: Design and Fabrication of Submerged Cylindrical Laminates. Part I and Part II, Int. J. Solids & Structures, 1999, p. 1248–1295. [9] Eshelby, J. D.: The determination of the elastic field of an ellip- soidal inclusion, and related problems. JMPS, Vol. 11, 1963, p. 376–396. [10] Kachanov, L. M.: Introduction to Continuum Damage Me- chanics. Dordrecht (Netherlands): Martinus Nijhoff Pub- lishers, 1986. [11] Moreau, J. J.: Some numerical methods in multibody dynam- ics: Application to granular materials. Eur. J. Mech. Solids Vol. 13, No. 4, 1994, p. 93–114. [12] Onck, P., van der Giessen, E.: Growth of an initially sharp crack by grain boundary cavitation. JMPS, Vol. 47, 1999, p. 99–139. [13] Procházka, P. P., Kugblenu, M.: Application of certain dis- crete element methods to the problems of fracture mechanics. Acta Polytechnica, Vol. 42, No. 4, p. 57–67. [14] Procházka, P. P., Šejnoha, M.: Development of debond region of lag model. Computers & Structures, Vol. 55, No. 2, 1995, p. 249–260. [15] Triantafillidis, N., Schraad, M. W.: Failure in aluminium honeycombs. JMPS, Vol. 47, 1999, p. 1093–1124. [16] Zienkiewicz, O. C.: The Finite Element Method. New York: McGraw-Hill Book Company, 1977. Appendix � � log d log d A B A B r x x a x b a b b a x x x x 1 1 2 2 1 2 21 2 2 2 � �� � � � � � �log arctan b a � � � (A.1) x x r x x x i j d A B 2 1� � 1. � � � � � � � � � � x x a x a b x a x a x x 1 2 1 2 2 1d A B A Barctan arctan � � � � � 2. � � � � �� ax x a x a x a x a x x 1 B A d A B 1 2 2 1 2 2 2 2 log 3. � � � � � � � � � � a x a x a b x a x a x x 2 1 2 2 1 2 d A B B Aarctan arctan � � � � � 2 (A.2) x r x x x i d A B 2 1� � 1. x x a x x a x a x x 1 B A d A B 1 2 2 1 2 2 2 2 � � � �� log 2. a x a x x a x a x x 1 2 2 1� � � � � � � � � � d A B B Aarctan arctan (A.3) © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 49 Acta Polytechnica Vol. 42 No. 4/2002 Fig. 7: Geometry for calculating boundary integrals d A B B Ax r a x a x a x x 1 2 1 � � � � � � � � � � � � arctan arctan (A.4) x x x n r x x x a r x x x x x i j k j i kd d A B A B 2 1 2 1� �� � see (A.2) Prof. Ing. RNDr. Petr P. Procházka, DrSc. phone: +420 224 354 480 fax: +420 224 310 775 e-mail: petrp@fsv.cvut.cz Ing. Michael Kugblenu Dept. of Structural Mechanics CTU in Prague, Faculty of Civil Engineering Thákurova 7 166 29 Prague 6, Czech Republic 50 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 42 No. 4/2002