Development of a Formalism of Movable Cellular Automaton Method for Numerical Modeling of Fracture of Heterogeneous Elastic-plastic Materials

A general approach to realization of models of elasticity, plasticity and fracture of heterogeneous materials within the framework of particle-based numerical methods is proposed in the paper. It is based on building many-body forces of particle interaction, which provide response of particle ensemble correctly conforming to the response (including elastic-plastic behavior and fracture) of simulated solids. Implementation of proposed approach within particle-based methods is demonstrated by the example of the movable cellular automaton (MCA) method, which integrates the possibilities of particle-based discrete element method (DEM) and cellular automaton methods. Emergent advantages of the developed approach to formulation of many-body interaction are discussed. Main of them are its applicability to various realizations of the concept of discrete elements and a possibility to realize various rheological models (including elastic-plastic or visco-elastic-plastic) and models of fracture to study deformation and fracture of solid-phase materials and media. Capabilities of particle-based modeling of heterogeneous solids are demonstrated by the problem of simulation of deformation and fracture of particle-reinforced metal-ceramic composites.


INTRODUCTION
modern theoretical approach used for description of heterogeneous solid is based on one of two concepts of representations of the medium: continuum or discrete.These concepts were suggested as memorandums to Paris Academy of Sciences by Cauchy and Navier, correspondingly, in the early 19th century.Notwithstanding the fact that discrete concept reflected discrete structure of matter more correctly at different spatial scales, during the following one and the half century the continuum concept prevailed.This was due to the fact that this concept gives the advantage of analytical description.It is necessary to note, that most of computational techniques used to model the behavior of condensed media is based on the continuum concept.The best known within this class of methods are finite element and finite-difference methods [1,2].At the present time meshless algorithms for numerical solution of equations of continuum (particle-in-cell method [3], SPH [4,5], SPAM [6,7]) are widely developed and applied as well.Numerical methods belonging to continuum concept has been proven quite efficient in solving problems of deformation of complex heterogeneous media of various nature.However, there exists a wide range of problems in which loading of a medium involves multiple fracture and slip of fragments, intense mass transfer, including effects of mass mixing, etc. Examples of these problems are flow of granular and loose media (powders, sands, soils), compaction of powder mixtures, processes in the contact spots of friction pairs, etc.These problems are too difficult to solve by methods of the continuum concept.
A A promising class of numerical methods that "genetically" fit for simulation of fracture of solids and mass mixing are particle-based methods.These methods are widely used to study mechanical or thermomechanical response of solid-phase and liquid-phase systems at various scales up to the atomic scale.Note that the term "particle methods" is currently understood as a collective term and is referred to quite diverse numerical techniques including methods based on the formalism of the discrete concept in mechanics of solids as well as meshless numerical methods of continuum concept.Moreover, nowadays some modern realizations of conventional numerical methods (such as particle-finite element method [8]) are also referred to as particle methods.The following consideration will concern "conventional" particle methods.
In the framework of "conventional" particle methods simulated material is considered as an ensemble of interacting particles (elements) having finite size.Evolution of an ensemble is defined by solution of the system of Newton-Euler motion equations: where i R  and i   are radius-vector and rotation angle of the particle i; m i and ˆi J are particle mass and moment of inertia; i F  is the total force applied to the particle i from the surroundings; ij n F  and ij t F  are forces of normal and tangential interaction of considered element i with neighbor j, ij M  is momentum of force, Ni is a number of neighbors.
Conventionally only nearest neighbors of particle i are taken into account in (1).It is seen from ( 1) that "macroscopic" (integral) properties of ensemble of particles are defined by the structure and parameters of potential (potential forces) of element interaction.
It should be noted that the particle methods have their origin in the well-known molecular dynamics methods (MDM) [9] applied to study of the response of a medium at the atomic scale.At the same time, capabilities of atomic description of the behavior of the solid body on spatial and temporal scales that are of interest for engineering applications are severely limited.This led to the development of numerical methods for meso-and macroscopic description of a medium (particle methods) on the base of MDM.In the framework of these methods structural elements (particles) have finite size and hence interaction with the immediate surroundings only is taken into account.
The best known representative of conventional particle-based methods is the discrete element method (DEM) [10,11].Its basic foundation was independently proposed by P.A. Cundall (distinct element method [10]) and Greenspan (particle modeling [12]) at the end of 1970s.At the present time these two similar methods have been developed into a quantity of methods and models, which are referred to collectively as DEM.
A distinctive feature of DEM in comparison with other particle-based methods is the presence of predefined initial shape of elements.Element shape can change (reversible or irreversibly) as a consequence of loading.This feature determined preferential development of the formalism of DEM towards correct modeling of flow of loose or weakly bonded porous materials [13][14][15][16][17].In the numerical description of such systems, an adequate accounting of geometric characteristics of discrete elements imitating particles of loose medium and peculiarities of contact interaction of elements is of crucial importance.In a general case particles of the simulated medium (discrete elements) are considered as super-quadrics [15] or polygons [18] rather than as disks or spheres (these shapes are used within the simplest models of loose material).Central forces ij n F  therewith can generate torque, which is impossible in the interaction of disks or spheres.At the same time, far less attention is paid to the analysis of structural form of potential (or potential force) of element interaction in loose medium.Indeed, the vast majority of models of interaction between discrete elements is based on the use of approximation of pair-wise elastic or elastic-plastic interaction.Such simplification is related, in particular, with the need to model long-time evolution of an ensemble of a large number of elements (up to tens millions particles in 3D problems).It should be noted that the range of problems that would be correctly solved with use of the approximation of pair-wise interaction between elements is limited by description of deconsolidated (including decompacted) materials and media with a large enough percentage of free volume (particulate media).In such systems the influence of lateral spread of deformed elements can be neglected.Therefore, the main areas of application of models with pair-wise interaction of discrete elements of complex shape are the theoretical study of flow (including transportation), mixing, segregation and compaction of granular materials in industrial and laboratory plants as well as under natural conditions (modeling of landslides and avalanches).
At the same time, the success of DEM in description of mechanical processes in consolidated solids (metals, low porous ceramic materials and rocks) is less visible.Hereinafter the term "consolidated solid" means low porosity or pore-free material retaining its shape under unconfined conditions.The characteristic features of the interaction of discrete elements, whose ensemble models consolidated material, are taking into account the resistance of element pairs both to compression and tension as well as limitation of the value of shear resistance ( ij t F ) by adhesion/cohesion strength rather than by coefficient of friction.Particular form of discrete elements modeling fragments of consolidated material or medium does not play a fundamental role.Conventionally the disk (in 2D problems) or spherical (in 3D problems) form of the element is used.At the same time the structural form of normal and tangential potential forces of discrete element interaction is of crucial importance for modeling consolidated materials.In particular, the use of conventional models of pair-wise potential interaction can lead to a series of artificial manifestations (effects) of response of the ensemble of elements that are not inherent to modeled medium.Most important of them are:  strongly pronounced dependence of macroscopic mechanical properties of ensemble of discrete elements on packing type (close, square, stochastic,…);  problems in realization of desired ratio between macroscopic elastic moduli (shear and bulk moduli, Young modulus and Poisson ratio and so on), especially in case of regular packing;  problems in correct simulation of irreversible strain accumulation in ductile materials, whose plasticity is provided by mechanisms of crystal lattice scale.At present there are several methods of partial solutions of these problems.They are associated with generation of stochastic dense packing of nonuniform-sized elements [19], definition of particle interaction constants using lattice approximation of continuum [20] and so on.However, the capabilities of these solutions are strongly limited, because they do not remove the main limitation of the pair-wise approximation in description of element interaction, namely the neglect of the change of discrete element volume under loading.This problem may be solved with use of many-body interaction forces.At present time different authors develop at least several approaches to description of the response of consolidated solids by ensemble of discrete elements with use of many-particle interaction models [20,21].However, the ability of these models at present is limited to description of elastic-brittle materials.Therefore the present paper is devoted to the development of a general approach to building many-body forces of discrete element interaction to simulate deformation and fracture of consolidated heterogeneous media with various rheological characteristics.The proposed approach is realized within the framework of movable cellular automaton method (MCA) [22][23][24][25].MCA method is a hybrid computational technique combining mathematical formalisms and capabilities of DEM [10,11,19] and conventional concept of cellular automata [26][27][28].Originally MCA method was developed to study complex coupled processes in solids including deformation and fracture, phase transitions, chemical reactions and so on [22,29].In this case, to describe the mechanical interaction of elementary structural units (cellular automata) the basic postulates and relations of DEM are used.Therefore in modeling "pure" mechanical problems MCA method can be considered as an implementation of DEM, which has the following principal difference from other implementations.During the construction of the mathematical formalism of MCA method general structural forms of the central and tangential potential forces of interaction of movable cellular automata were derived.Central interaction between movable automata has many-body form, whereas tangential one has pair-wise form [24]. Thus, the interaction between movable cellular automata originally supposed to be many-particle.Note that conventional formalism of discrete elements does not have any supposition about the form of interaction potential/force.In this paper consideration will be conducted by the example of two-dimensional problem statement, which provides the possibility of movement of three-dimensional objects in one plane only.Construction of many-particle interaction in the three-dimensional problem formulation is similar.Note also that the following approach and models are applicable not only for MCA method, but for a wide range of implementations of discrete element method.Therefore both terms (movable cellular automata and discrete elements) will be used in the text.

MAIN PARAMETERS OF INTERACTION OF DISCRETE ELEMENTS.
n the framework of proposed approach to modeling consolidated solids the description of the interaction of discrete elements (or movable cellular automata) is based on the use of two types of states of a pair of interacting elements.They are associated with presence (linked state) and absence (unlinked state) of chemical bond (or I cohesion/adhesion) between elements (Fig. 1).The feature of the central interaction of linked elements is the presence of resistance to both compression and tension.At the same time, the interaction of unlinked elements relates to the contact type and involves only the resistance to compression.Another difference of element interaction in linked and unlinked pairs is the fact that the magnitude of the tangential potential force is limited by the bonding strength in linked state, and by the dry friction force for contacting (unlinked) pairs.Thus, the ensemble of linked discrete elements simulates consolidated solid.In this case, discontinuities in such solid (damages, cracks, pores and so on) can be modeled by specifying unlinked pairs or by removing discrete elements.Switching between the linked and unlinked states of the pair is based on chosen criteria.Specific form of the criteria of direct (linkedunlinked) and reverse switch is defined by the nature of the materials modeled by discrete elements.In general it can be noted that switching criteria are functions of spatial and force parameters of the mechanical interaction between the elements.More details on this aspect of the interaction between discrete elements are discussed in Section 4. In framework of the conventional models (1) the mechanical interaction between particles is divided into normal or central (along the line connecting the centers of mass of elements) and tangential constituents.Each of these components is controlled by the corresponding spatial parameters.Normal interaction is determined by the value of element-element overlap h ij : where r ij is the current distance between the centers of mass of the elements i and j, 0 ij r is the original distance (in undeformed state), di and dj are the sizes of elements.Linked pair of discrete elements can resists both to compression and tension, therefore the value of h ij could be positive (tension of the pair) or negative (compression).In general, the discrete elements i and j can be characterized by different material properties, so that the contribution of each of them to hij may be different.In this regard, the concept of the distance between the center of mass of the element and the central point of plane of interaction (or, using the traditional terminology, the contact point) is introduced to the model: where q ij and q ji are the corresponding distances (Fig. 2).Initial values of q ij and q ji are d i /2and d j /2 correspondingly.
Figure 2: Parameters of spatial relation of the pair of discrete elements i and j: distance between mass centers (r ij ) and distances from mass centers of interacting elements to the center of plane of interaction (q ij and q ji ).
For convenience hereinafter spatial parameters of central interaction of discrete elements will be considered in reduced units (deformations): where symbol  hereinafter indicates increment of corresponding parameter during one time step t, ij is central strain of the pair i-j, variables  i(j) and  j(i) are central strains of discrete elements i and j in the pair (in the general case  i(j)  j(i) ).Tangential interaction is determined by the relative shear displacement shear ij l of the elements in the pair.The value shear ij l is calculated in incremental fashion taking into account the rotation of elements [19,22]: where tan g ij V is the tangential component of the relative velocity vector ij V  for the elements i and j ( ij and j V  are velocity vectors of the centers of mass of the elements), i and j are angular velocities of the elements (in considered two-dimensional approximation they actually are signed scalars).Note that the accounting of angular velocities in ( 5) is required for tracking the rotation of the plane of interaction [19,22].As in the case of the central interaction, the contribution of the elements i and j in the shear deformation of the pair i-j can be different.By analogy with the central strains (4) the shear angles of the discrete elements i and j in the pair i-j are introduced: where ij is shear angle for the pair i-j, variables i(j) and j(i) are shear angles of discrete elements i and j in the pair (in the general case  i(j)  j(i) ).Rotation of discrete elements can lead to "bending" of linked interacting pairs (Fig. 3).This type of relative motion of surfaces of elements i and j is not taken into account by expression (5).Nevertheless such type of relative motion of "anchor surfaces" must be accompanied by appearance of special moment of resistance force.So, the special torque gear K  directed against pair bending (Fig. 3) and having opposite signs for two elements of the pair is introduced: Value of the torque gear  is defined by bending angle gear ij  , which is calculated by analogy with shear ij  : By analogy with the central (4) and shear (6) strains the bending angles of the discrete elements i and j in the pair i-j are introduced: where Note that in the "conventional" models of the interaction of discrete elements do not take into account "bending" of the pair of discrete elements.This is perfectly valid in describing the contact interaction of unlinked elements.However, when considering linked pairs, accounting of this effect may be crucial.It follows from ( 5), ( 6) and (9) that strains of elements i and j in the pair i-j differ from each other.The rule of strain distribution in the pair is inseparably linked with the expression for element interaction forces and will be discussed in the following Section.By analogy with strains the forces of central ( ij n F ) and tangential ( ij F  ) interaction of discrete elements i and j will be considered in reduced (specific) units: where  ij and  ij are normal and shear stresses (specific forces) correspondingly, S ij is an area of the plane of interaction of elements (contact area).
Note that the definition of Sij is different for the description of consolidated and granular materials by ensemble of discrete elements.Thus, in the case of granular material the initial value of S ij in the undeformed state ( 0 ij S ) is defined by shape of the elements i and j and by their mutual orientation.For example, 0 0 ij S  for the elements having a disk (2D case) or spherical (3D case) shape.With increase of compression the value of contact area grows up as a function of hij: This kind of dependence is determined by the geometrical characteristics of both elements and their mutual orientation.At the same time, under the simulation of consolidated material the parameter 0 ij S is obviously finite and determined by the ratio of sizes of the elements i and j and by features of the local packing of discrete elements in the neighborhood of the pair i-j.As a rule, the 0 ij S is defined to minimize the volume of voids in simulated solid (note that in bonded-particle model [19] this is realized through introduction of the special interfacial zone between two bonded (linked) elements and by assigning its geometrical characteristics).Peculiarities of the procedure of assigning the initial area of plane of element interaction lead to the fact that in most models of consolidated solids a simplified form of discrete elements is used (disks or spheres).In two-dimensional approximation discrete elements are considered as disks with diameters di in the plane of motion and the same value of height h for all of them.Note that in the case of a regular packing of elements with the same size d i =d the value 0 ij S is determined by type of packing.For example, for a close packing 0  3 ij S dh  , for a square package 0 ij S dh  (Fig. 4).Change of S ij during deformation of the pair of discrete elements i and j is determined by approximations of applied model of element interaction.In particular, at small values of pair strain .However, at larger pair strains the value of S ij has to be considered as a function of the deformation process in the pair i-j.

GENERAL FORMALISM OF PROPOSED APPROACH TO BUILDING MANY-PARTICLE INTERACTION OF DISCRETE ELEMENTS
uthors propose a general approach to building many-body forces of discrete element interaction to simulate deformation and fracture of consolidated heterogeneous materials.The structural form of these forces is borrowed from the form of interatomic forces calculated on the basis of embedded-atom method.In the framework of embedded-atom model [30] the general expression for potential energy of atom i contains a pair interaction potential  as a function of distance r ij between atoms i and j and a "density-dependent" embedding function F (here it depends on electron charge density i where is a sum of contributions of neighbors j to local value of density at the location of atom i.
By analogy with this expression the following general form of notation of the expression for the force i F  acting on discrete element i from surroundings is proposed [23,31]: This force is written as a superposition of pair-wise constituents ij pair F  depending on spatial position/displacement of element i with respect to nearest neighbor j and of volume-dependent constituent i F   connected with combined influence of nearest surroundings of the element.When simulating locally isotropic materials/media with various rheologies the volume-dependent contribution i F   can be expressed in terms of pressure P i in the volume of discrete element i as follows [23]: where Sij is square of area of interaction (contact) of elements i and j, ij n  is a unit vector directed along the line between mass centres of considered elements, Ai is a material parameter.
In such a formulation the right part of the expression (12) can be reduced to the sum of forces of interaction in pairs of elements and divided into central where pair ij  and pair ij  are specific values of pair-wise components of interaction force.In fact ij and ij have the meaning of specific forces of response of the element/automaton i to the impact of the neighbour j.Taking into account the need to implement Newton's third law for interacting pairs of discrete elements (ij=ji and ij=ji) the expression (15) has to satisfy the following equality: The equality (15) provides the basis for calculation of contributions of elements i and j ( i(j) and  j(i) ,  i(j) and  j(i) ) in the total pair strains  ij and  ij .Note that although the right part of the expression ( 14) formally confirms to notation of element interaction in conventional DE models (1) [16,19,32], their fundamental distinction consists in many-body form of central interaction of discrete elements in the proposed model.It is seen from ( 14)-( 15) that an important problem in building many-particle interaction is definition of local value of pressure (Pi) in the volume of discrete element.Authors propose to use an approach to calculation of pressure Pi (or, what is the same -of mean stress) in the volume of the element i that is based on the computation of components of average stress tensor in the volume of the element [19,33].In terms of central ( ij n F ) and tangential ( ij t F ) interaction forces the component i   of average stress tensor in the volume of element i can be written as follows [19, 24, 25]: where , = x,y,z (XYZ is a laboratory system of coordinates),  i is a current value of the volume of element i,   are projections of unit-normal and unit-tangential vectors onto the X-axis of lab coordinates.
In considered two-dimensional problem statement this expression can be rewritten as follows: cos cos cos sin where , = x,y (XY is a plane of motion), ij,  is an angle between the line connecting mass centres of interacting elements i and j and axis  of laboratory system of coordinates (Fig. 5  Calculated in this way the stress tensor components can be used to determine the pressure in the volume of discrete element: Note that calculated values of average stress tensor components can be used to determine other tensor invariants as well, for example stress intensity: It follows from ( 1), ( 14)-( 16) that the central problem in the framework of proposed approach to building many-body interaction of discrete elements (or movable cellular automata) is to determine expressions for ij and ij, which provide necessary rheological characteristics of mechanical response of ensemble of elements.Analysis of relationships ( 1), ( 14)-( 16) leads to the conclusion that expressions for interaction forces could be directly reformulated from constitutive equations of considered medium (equations of state).
Below is a derivation of such expressions for locally isotropic elastic-plastic materials.

DISCRETE ELEMENT INTERACTION FOR MODELING CONSOLIDATED ELASTIC-PLASTIC MEDIUM
Linearly elastic medium tress-strain state of isotropic linearly elastic medium is described on the basis of generalized Hooke's law.The following notation of this law will be used in the paper: where , = x,y,z;   and   are diagonal components of stress and strain tensors;   and   are off-diagonal components; It can be seen that the form and the matter of expressions (19) for diagonal and off-diagonal stress tensor components are analogous to expressions (15) describing normal and tangential interaction of discrete elements.This leads to the simple idea to write down expressions for force response of automaton i to the impact of the neighbor j by means of direct reformulation of Hooke's law relationships: where Gi and Ki are shear and bulk elastic moduli of material filling the element i, i(j) and i(j) are central strain and shear angle of element i in the pair i-j (see Section 2), mean stress i mean  is calculated using (17).
It is necessary to note that double shear modulus (2G i ) is used in the second expression of ( 20) instead of G in (19).This feature is concerned with the fact that relative tangential displacement of discrete elements leads to their rotation.Initiated rotation of the elements decreases twice the value of relative shear displacement in interacting pairs.So, to match "macroscopic" (integral) shear modulus of ensemble of discrete elements simulating consolidated solid with corresponding shear modulus G double value of shear modulus (2G i ) is used in the second expression of (20).Proposed relationships (20) for forces of element response to the impact of the neighbor j are not arbitrary.It is easy to verify that substituting proposed expressions (20) for respond force in ( 16) automatically provides implementation of the Hooke's law for components of average stress ( i   ) tensor in the volume of element i. Detail description is presented in Appendix A. Proposed relationships (20) make it possible to calculate central and tangential interaction of discrete elements, whose ensemble simulates isotropic elastic medium.Taking into account the need to implement Newton's third law for S interacting pairs of discrete elements (15') and the need to distribute relative displacement of elements in the pair the expressions for specific interaction forces can be written as follows: and Here, relations for calculating the central and tangential interaction forces are written in incremental fashion (in hypoelastic form); "cur" and "pre" upper indexes mark values of specific reaction forces at the current step of integrating the equations of motion of discrete elements (or cellular automata); mean stress increments i mean   and j mean   are taken from previous time step or determined with use of predictor-corrector modification of a numerical scheme.Equations ( 21)-( 22) are first solved for strain increments . Found values of strain increments are then substituted in ( 21)-( 22) to calculate current values of specific forces cur ij  and cur ij  .
The above relations have a general (three-dimensional) form.In considered two-dimensional problem statement the approximations of plane stress (zz=0) or plane strain (zz=0) state are used.In the model under consideration, this means that the average stress tensor components i xz where  i is a Poisson ratio.
Resistance of the pair i-j to bending (Fig. 3) is not included in the expressions ( 21)- (22).So, the form of the dependence has to be assigned extra.In the present model this is done by analogy with (22): Testing of the proposed model of element interaction by the example of elastic wave propagation through ensemble of linked elements showed that use of shear modulus G as an elastic modulus of resistance to bending in (24) provides the best agreement of simulation results with analytical solutions and results of simulation by finite difference method.Relations ( 21)-( 24) make it possible to simulate fragments of a locally isotropic linear elastic medium by an ensemble of linked discrete elements (or movable cellular automata).

Elastic-plastic medium
An important advantage of the proposed approach to building many-body interaction of discrete elements is a capability to realize various models of elasticity and plasticity within the framework of DEM/MCA.In particular, a model of plastic flow (incremental plasticity) with von Mises criterion of plasticity was implemented to simulate deformation of locally isotropic elastic-plastic medium.
For this purpose, radial return algorithm of Wilkins [2] was adopted to particle-based approach.Conventionally, this algorithm is formulated in terms of the stress deviator D (Fig. 6): ˆD D M M is a coefficient of stress drop (stress scaling), D is a stress deviator after solution of elastic problem at the current time step, D  is a scaled stress deviator.
Being written in terms of stress, for components of average stress tensor in the volume of discrete element i the algorithm of Wilkins can be presented in the following form: where  are scaled (returned) components of average stress tensor; i   are stress tensor components, which result from solution of elastic problem ( 21)-( 23) at the current time step;  is current radius of von Mises yield surface for the element i; int i  is calculated with use of expression (18) after solving elastic problem at the current time step;   is the Kronecker delta.The main problem in realization of the algorithm of Wilkins within the framework of DEM/MCA is formulation of correcting relations for element interaction forces that provide implementation of necessary conditions of the algorithm [2].By analogy with the elastic problem the expressions for scaling specific central and tangential forces of response of the element i to the impact of the neighbor j were derived by direct reformulation of relations (26) for average stress: where ij   and ij   are scaled values of specific reaction forces.
As is shown in Appendix B, substitution of (27) in expression (16) for average stress tensor automatically provides reduction of its components to yield surface for the element i.This gives possibility to correctly simulate plastic deformation in the volume of discrete elements.Note that independent use of the expressions (27) for interacting elements i and j can lead to unequal values of respond ) in the pair i-j.In view of the need for implementation of Newton's third law, the current values of element interaction forces in (10) are calculated on the basis of the following proportion: In considered two-dimensional problem statement scaling of specific forces  ij and  ij and stress i zz  has peculiarities for approximations of plane strain and plane stress state.In the first case these variables are scaled using the expressions (26) for i zz  and ( 27) for specific forces.In the second case the special iterative procedure proposed by Wilkins [34] is adopted to particle-based approach (see Appendix C).By analogy with the case of elastic problem, law of scaling of resistance of the pair i-j to bending has to be defined extra.In the present model this was done by analogy with scaling of shear resistance force  ij : Rheological properties of material of discrete element i are defined through assigning constitutive relation (when applied to movable cellular automaton, it is called as "mechanical response function of automaton/element" [22][23][24][25]).Current value of int i  could be calculated incrementally using known values of int i  after solution of elastic problem at the considered time step (n+1) and at the end the previous time step n: where     remains unchanged throughout the stress return procedure (this is applied both to central and shear strains of discrete elements i and j in the pair i-j).

Calculation of current values of element volume and square of area of interaction of the pair
In addition to forces  ij , and  ij , acting on the surface of discrete element i, important constituents of the expression (16) for the components of average stress tensor are element volume i and squares of areas of interaction of the element with neighbours Sij.The current values of these variables can be found using the average strain tensor components i   in the volume of the discrete element i.The values of components i   can be found directly in terms of pair strains  i(j) and  i(j) of element (see Appendix A) or through the components of the average stress tensor i   .In the last case, the specific relations between i   with i   are determined by considered rheological model of the medium.In the framework of above described two-dimensional model of elastic-plastic interaction of discrete elements with von Mises criterion of plasticity these relations should be presented in hypoelastic form: are stress increments calculated after solving the elastic task at the current time step.


 are used to define current volume of the discrete element i: where 0 i  is "initial" volume of the element (in undeformed state).Determination of current value of S ij is more complicated problem due to possible significant element shape distortion under large loads.The following approximation to estimate the value of Sij is suggested.It is based on definition of local values of strain tensor components at the area of interaction (contact area) of considered pair of discrete elements i and j (hereinafter denote this tensor as ij   ) by linear interpolation of corresponding components of average strain tensors ( i   and j   ) to the central point of this area: Components of ij   tensor are transformed from the laboratory system of coordinates to the instantaneous local coordinate system XY of the considered pair (Fig. 7)    and ij zz  thus defined in local coordinate system are used for calculation of the current value of square of the area of pair interaction: Here 0 ij S is the initial value of square corresponding to the pair of undeformed elements i and j.

Organization of the step of numerical integration of motion equations
Thus, numerical algorithm realizing the proposed model of elastic-plastic interaction of discrete element (or movable cellular automaton) i with neighbours j includes the following main stages: i. Calculation of ij, ij and gear ij K at the current time step (n+1) in elastic approximation according to Eqs. ( 21)- (24). with use of expressions ( 16)-( 18), ( 23), (30).iii.Examination of exceeding the yield surface by element i stress state controlling parameter int i  .Calculation of the coefficient M i if necessary.iv.Carrying out the procedure of stress returning to yield surface for element i with use of expressions ( 26)-( 27), (29)

Calculation of corresponding values of
Calculation of scaled values of average stress tensor components and invariants.
v. Calculation of forces and torques of element interaction.vi.Calculation of new values of automaton volume i and squares of areas of pair interaction Sij with use of expressions ( 32)-( 34) vii.Integration of Newton-Euler equations of motion of discrete elements (or movable cellular automata).The equations of motion (1) of discrete elements are numerically integrated with the use of the explicit integration scheme (for example, with use of well-known velocity Verlet algorithm) [35] modified by introducing a predictor cycles for estimation of i   at the current step.The need of the numerical scheme modification is caused by the fact that specific forces of the central interaction of discrete elements  ij at the current time step are computed using mean stress i mean  (21).At the same time these forces are used to define mean stress according to (16).To solve this problem of numerical integration the predictor cycles are used for the calculation of the elastic task.At the initial predictor cycle increment of mean stress is calculated from the previous time step:

MODELING FRACTURE AND "HEALING" WITH PARTICLE-BASED FORMALISM
ne of the main advantages of particle-based methods in mechanics is the feasibility of direct simulation of fracture (including multiple fracture) of material and bonding of fragments through changing the state of a pair of particles ("linked" pair  "unlinked" pair, Fig. 1).The criterion for pair state switching is normally the ultimate value of interaction force or the ultimate value of relative displacement [19].The developed approach to the description of interaction of movable cellular automata (or discrete elements) in the many-particle approximation makes it possible to apply various multiparametric "force" fracture criteria (Huber-Mises, Drucker-Prager, etc) as element-element bond fracture criteria.

The model to calculate debonding of the pair (linkedunlinked transition)
It's well known that crack formation is a fundamentally brittle and extremely localized phenomenon.Fracture is connected with rupture of interatomic bonds and spatial diversity of atomic layers.That is realized under the influence of local, rather than the average stress.Thus, the used criterion of pair bond breakage between linked elements i and j must be expressed in terms of local variables of the interaction of elements, in particular, through the forces  ij and  ij acting in a pair.In this regard, the following approach to implement parametric fracture criteria as criteria of pair bond breaking (linkedunlinked transition of the state of the pair) is suggested.In the framework of classical formalism of discrete elements pair bond breakage occurs on the surface of their interface (at the area of interaction of the pair, in other words, O at the contact area).Therefore applied "force" failure criterion (for example, the criterion of Drucker-Prager) must be calculated at the area of interaction of elements using the local stress tensor components identified at this area.By analogy with strain tensor ij   (see Section Calculation of current values of element volume and square of area of interaction of the pair) this local tensor will be denoted as ij   .In the local coordinate system XY of the pair i-j (Fig. 7) components ij y y    and ij x y    of this stress tensor are numerically equal to specific forces of central ( ij ) and tangential ( ij ) interaction of the elements (these forces are applied to the contact area S ij ):    ) of the local stress tensor can be determined on the basis of linear interpolation of corresponding values for elements i and j ( i ) to the area of interaction: where i      and j      are components of average stress tensor in the volume of elements i and j in the local coordinate system XY of the pair (these stresses result from transformation of average stresses i   and j   to local coordinates).
Components ij      , thus defined, can be used to calculate necessary invariants of stress tensor which then can be used to calculate current value of applied criterion of pair fracture.Below the examples of bond breaking conditions for the pair i-j with use of Huber-Mises-Hencky and Drucker-Prager criteria are shown: where  c is the corresponding threshold value for considered pair (value characterizing strength of cohesion/adhesion), When using the explicit scheme of integration of motion equations (1) the value of time step t is limited above by a quantity associated with the time of propagation of the sound through the volume of element: where V elast is velocity of the longitudinal elastic wave in the material, 0<<1 -dimensionless coefficient (normally =1/51/4).In conventional DEM-based models breakage of bond in pairs of discrete elements is carried out for one time step t.It this case the local value of the crack growth rate (defined as the ratio of the linear size of the area of pair interaction to t) is higher than sound velocity in the material.Obviously, such a model of bond breakage is an idealized because virtually suggests that the spatial separation of atomic layers occurs uniformly over the whole surface of interaction of elements.In modern models of fracture mechanics it is supposed that fracture has a discrete character, that is, the crack length increases by "quanta".In different models these quanta of length are called as "fracture quantum" [36,37] or "process zone" [38,39].Thus, in conventional DEM-based models of fracture it is implicitly assumed that the linear dimension of the surface of interaction of elements is equal to the step of crack growth.At the same time, the value of "fracture quantum" is a material parameter and its scale normally is nanometers or fractions of nanometer [39,40].In the most tasks the size of discrete elements is significantly larger than this value.Therefore, in most cases, the model of bond breaking in pair of discrete elements within one integration step t means the overestimation of the rate of physical process of fracture of the material modeled by the pair of linked discrete elements.
In accordance with the foregoing, the following approach to a more accurate description of the dynamics of crack growth is suggested.In this approach, it is assumed that breaking of the bond (linkedunlinked transition of the state of the pair) is a time-space distributed process.This process is technically expressed through change of the dimensionless coefficient . This coefficient has the meaning of the portion of linked part of the contact area Sij.In this case the square of linked part of the contact area in the pair i-j is the square of unlinked part of this area.Thus, in this approach, the dynamics of bond breaking in a pair of elements is expressed by the dependence   ij link k t , where ij link k decreases from initial value 1 (totally linked pair) to final value 0 (totally unlinked pair).Depending on the size of the discrete elements and features of the internal structure of fragment of the material, which is simulated by the discrete element (in particular, the presence of pores, damages, block structure) the stable or unstable crack growing model can be applied to describe the breakage of bond in the pair.In the first model, the process of fracture develops in accordance with predefined dependence , where int ij  is a pair strain intensity (value int ij  can be calculated, for example, using the components of the tensor ij   defined by analogy with ij   ).In the simplest case this dependence can be considered as a constant:  .Note also that the additional conditions for decrease of ij link k (i.e.conditions of crack advance through the area of pair interaction) at the current time step are: i. Exceeding the fracture criterion threshold at this step.ii.int 0 dk d has to be assigned for each pair of materials filling linked elements.In the second model (the approximation of unstable crack growth) it is assumed that if the fracture criterion threshold is exceeded, a crack begins to grow spontaneously according to some specified law , where t is time.
In the simplest case this dependence can be considered as a constant ( 0 ij link dk dt const   ), that means that crack advances through pair interaction area with constant velocity V crack .The value of V crack is a predefined (model) parameter, which reflects the features of rheology of the interface between the materials filling interacting elements.In particular, for brittle materials V crack can be set to be equal to transverse sound velocity, while for ductile materials its value, obviously, should be significantly smaller.Note that the model of unstable crack growth has to be used together with requirement of exceeding the fracture criterion threshold at each time step during crack advance through the area of pair interaction.Model parameter V crack has to be assigned for each pair of materials filling linked elements.

Possibility of unlinked to linked transition
Chemical bonding of contacting unlinked automata is imitated by unlinkedlinked transition of the pair state (this transition is interpreted as formation of new chemical bond or recovery of previously broken one).Note that such transition describes "healing" of partially linked pairs of elements ( 0 1 ij link k   ) as well.Physical mechanisms of "integration" (linking) of independent material fragments into a consolidated piece could be different and include cohesion (or adhesion) of smooth and pure enough contacting surfaces under compression, "welding" of contacting surfaces under the condition of compression and intensive friction, interpenetration ("mixing") of surface layers of contacting material fragments as a result of strong compression and intensive shear (torsion) under constraining side boundary conditions, healing of nano-and microscopic damages and cracks and so on.Therefore specific form of bonding/linking criterion is defined by physical peculiarities of contact interaction of automata as well as by chemical composition and structural features of interacting surfaces.Due to a surface roughness the process of "integration" of surface layers of contacting automata is gradual.The degree of "integration" of surface layers increases with increasing the values of variables of loading (normal load, shear load, plastic work of deformation).To take into account these features, the following model of linking/bonding of unlinked and contacting automata is suggested.The change of degree of "integration" of surface layers of totally or partially unlinked elements i and j is taken into account by means of change of the value of previously introduced coefficient ij link k (see Section The model to calculate debonding of the pair (linkedunlinked transition)).In the framework of this model the criterion of pair linking takes the form of dependence of ij link k on the variables of loading.The following simple and physically based criteria are proposed: 1. Linear dependence of ij link k on specific normal force ij: where  ij is a current value of specific central force in the pair i-j;  ij is an increment of specific central force during one time step t; min is the minimum (threshold) value of normal pressure in the pair (negative specific normal force) for linking; max (max>min) is another parameter having the dimension of specific force and regulating the slope of the dependence  is a maximum value of specific normal force previously achieved from the moment of linking beginning.
2. Linear dependence of ij link k on plastic work of deformation W ij taking into account specific normal force  ij : where W  is normalizing parameter depending on specific normal force ij: Here Wij is an increment of plastic work of deformation of the pair i-j during one time step t;  39)-( 40) are assigned (user-defined) values, which have to be determined for each pair of materials filling elements i and j.One of the problems with the use of the criterion (39) is the calculation of Wij.As a first approximation it can be considered as the sum of increments of plastic work of deformation of both interacting elements/automata: A specific expression to calculate increment of plastic work of deformation of the element is determined by the applied model of plasticity.In the case of the above model of plasticity with von Mises criterion the value of Wi can be calculated as follows: where total i A  and elast i A  are increments of total work of deformation and its elastic part correspondingly, n is a number of time step.In a general case criterion of pair linking has more complex form as it takes into account combined influence of central force, plastic work and other parameters of pair interaction including strain rate and/or relaxation times.Note that during time-distributed process of pair linking the condition of bond breaking in a pair can be achieved.Therefore, under certain conditions, local processes of linking and unlinking can proceed in parallel (unlinkedlinked transition).

MCA-BASED STUDY OF FEATURES OF DEFORMATION AND FRACTURE OF METAL-CERAMIC COMPOSITES
apabilities and advantages of the developed approach to the construction of many-body potentials/forces of interaction between particles (discrete elements or movable cellular automata) and implemented rheological models make it possible to study the response (including fracture) of heterogeneous materials and media of different nature.An important area of application of particle-based numerical modeling is investigation of the influence of the features of internal structure on the mechanical properties and fracture pattern in composite materials.This will be demonstrated by the example of MCA-based modeling of particle-reinforced metal-ceramic composite.Metal-ceramic composites are advanced representatives of the class of dispersion-reinforced materials, which have enhanced values of mechanical and service characteristics, such as strength, stiffness-to-weight ratio, crack growth resistance, wear resistance, fracture energy, ratio of thermal conductivity to thermal expansion coefficient, thermal stability and so on.This makes them very attractive for a wide application in various industries as materials for extreme operating conditions [41][42][43].At the present time working parts of machines and mechanisms operating under conditions of shock loading, abrasion, high temperatures and corrosive environments are made mostly of metal-ceramic composites on the basis of very hard and refractory compounds (carbides, nitrides, carbonitrides) with metallic binder (nickel and iron alloys are mainly used) [41,42,[44][45][46].This class of materials is fabricated from powder mixtures of the compounds using powder metallurgy methods [41,42,47].Mechanical and physical properties of the sintered metal-ceramic composite are determined, in addition to phase composition, by a number of structural factors.Some of them are conventionally observed: volume fraction, dispersion (including the size distribution), geometry and faultiness of reinforcing particles, structural-phase state of the metallic binder et al. [41,42,48,49].It should be noted that one of the key elements of the internal structure of the metal-ceramic composites are the interfaces between particles of refractory compounds and a metallic binder.The change in technological peculiarities of metal-ceramic composite fabrication (in particular, applying additional heat treatment of the composite) can vary the geometry (width) of interfaces as well as their structural and, consequently, the mechanical properties [50].The analytical models and numerical simulation are of great importance for studying the mechanical properties of composite materials and their connection with geometric and mechanical properties of inclusion/matrix interfaces.As a rule, the analytical models to study the effect of the properties of interfaces on the macroscopical mechanical properties of dispersion-reinforced composites are valid for elastic or viscoelastic approximations and a regular packing of very hard inclusions having a simple geometry.In this connection the predominant attention in solving such problems is given to the numerical simulation of deformation and fracture of composites taking into account their mesoscopical internal structure.Finite element and finite difference methods are most extensively used numerical techniques for this purpose [48,[51][52][53].With the application of these numerical methods the influence of various structural factors (geometry, volume fraction, size and spatial distribution of reinforcing particles) on the elastic and strength properties, deformation capacity and fracture toughness of dispersion-reinforced composites was investigated [48,52,54].At the same time, the influence of strength properties and the width of interphase boundaries interfaces on the mechanical characteristics of metal-ceramic composites is still under discussion.Thus in this paper, particle-based MCA method has been applied to study the influence of these factors on strength, ultimate strain and fracture energy of the composite under dynamic loading.A TiCparticle-reinforced Ni-Cr matrix composite (50 vol.%TiC) has been chosen as a model system classified as the metalceramic composite materials in which the particles are much stiffer than the binder.

Mesoscopical model of metal-ceramic composite
The TiC-particle-reinforced Ni-Cr matrix composite (50 vol.%TiC) was used as a prototype to create a numerical model of metal-ceramic sintered composite.A typical microstructure of Ni-Cr/TiC composite is shown in Fig. 8a.As can be seen from Fig. 8a, carbide component consists of particles with quite complex geometrical shapes.The mean size of reinforcing particles D mean is 2.71.2m, maximal particle size D max is 7.6 m, Fig. 8b.It is necessary to note that thermally activated diffusion processes during the course of sintering result in the formation of a transition zone (the region of variable chemical composition) at the ceramic/metal interfaces.The characteristic width (H if ) of the transition zone in considered composite is a few micrometers.Note that local profile of distribution of chemical elements is defined by the features of a local microstructure (including the distance to neighbouring TiC particle, shape of particle surface and so on).Thus, the metal-ceramic composite contains three main structural components: a metallic binder (Ni-Cr), integrated high-strength brittle inclusions (TiC) of mesoscopical scale (1-10 m in diameter) and transition zone "particle-binder" (the region of variable chemical composition), whose width can reach several micrometers.Mentioned above main features of the mesoscopical structure in metal-ceramic composite formed the basis for a twodimensional structural model developed in the framework of the numerical MCA method.Each of constituents is modelled by ensemble of movable cellular automata (discrete elements) with appropriate rheological parameters (thus the cellular automaton simulates a domain/fragment of an inclusion, a binder, or an interface zone).The size d of a cellular automaton is an assigned parameter of the model.In the presented model a requirement for the value of d was accepted to be a few times higher than the characteristic size of grains of the binder.In that case mathematical models of elasticity and plasticity of isotropic media can be correctly applied to describe the deformation of the material fragment with characteristic size corresponding to the automaton size d.The maximal value of the dimensional parameter d of the model determines a specification of TiC particle shape in detail."Optimum" definition range of the value d (the size of the discrete element) is 0.1-0.3m for the considered metal-ceramic composite (typical grain size of the metallic matrix in this composite does not exceed 0.1 m).
As an example of MCA-based structural model, Fig. 9 shows the structure of an idealized metal-ceramic composite with TiC particles having a nearly spherical shape and the average size D TiC of 3 m (the size of movable cellular automaton in this case is d=0.3 m).In this example the approximation of "monosize" distribution of TiC particles is used (sizes of the inclusions are uniformly distributed about the mean value D TiC ; the amplitude of deviation from the mean value of the particle size is 0.1DTiC).The volume fraction of carbide particles in the model composite is 45%, which is close to the corresponding value (50%) in the real composite.The spatial distribution of TiC inclusions in the shown idealized composite is slightly inhomogeneous (i.e.no pronounced clustering of particles), the characteristic distance between the surfaces of neighboring particles is about 1.5 micrometers that is 50% of DTiC.Depending on features of internal structure of real composite and on the level of detail of the model there are two ways to account spatial (including width) and physical-mechanical characteristics of the transition layer between reinforcing particles and a metallic binder: 1.The model of a "narrow" transition zone.In this model, it is assumed that the width of the interface is much smaller than the size of movable cellular automaton d.Capabilities of the model are limited by specifying certain values of strength parameters (adhesive strength) in pairs of movable cellular automata, one of which simulates a segment of the inclusion surface and the second one represents an adjoining (adjacent) segment of the binder.Such approach does not take into account geometric characteristics of the transition zone, the presence of a concentration gradient of the chemical elements, and hence the existence of a gradient in mechanical properties.The defining characteristics of a "narrow" interphase boundary are its strength properties.An assigned value of adhesive strength for the pair of dissimilar cellular automata corresponds to the effective strength of interphase boundary.Note that described model does not account for the features of the rheology of the transition layer and can be used effectively in the case where the width of this layer does not exceed the assigned size of movable cellular automaton d. 2. The model of a "wide" transition zone.In this model, it is assumed that the width of the interface is comparable or greater than the size of the movable cellular automaton d.Here particle/binder interface is regarded as an area of variable composition of chemical elements (Ti, Ni, Cr, C) and modeled by several layers of cellular automata.Physical and mechanical characteristics of these "transition" automata vary with a distance from particle surface into the volume of binder for a given law (in accordance with local chemical composition).Relations between a local content of chemical elements and the corresponding rheological properties of "transition" cellular automata must be determined on the basis of experimental data.The proposed model can be efficiently used in numerical simulation of composites in the case of the width of the transition zone higher than the assigned size d of a cellular automaton.To implement the explicit way to account the transition zone between particles and a binder in the numerical model of composite the special algorithm was developed.This algorithm uses information about a change in characteristic particle size during composite fabrication and resulting profile of distribution of chemical elements (for example, Ti) in formed transition layer.The output of the algorithm is a set of data about physical (density, concentration of chemical elements) and mechanical characteristics of movable cellular automata belonging to transition zone.These properties correspond to material at an appropriate distance from the surface of the inclusion.Example of a model metal-ceramic composite with explicit consideration of transition zone between particles (TiC) and the binder (Ni-Cr) is shown in Fig. 10.The main problem when using an explicit way of modeling the transition zone at particle/binder interface is the determination of rheological parameters and strength of transition areas (layers) that are at different distances from the particle surface.These parameters are difficult to determine in experiment.Therefore in the present study a simple (phase mixture) model was used to determine the mechanical properties of cellular automata that simulate transition zone.It is based on the assumption that the transition zone has a structure similar to the solid solution of titanium and carbon in the metallic binder.Their content in the binder decreases with increasing distance from the particle surface.In this model the rheological and strength parameters of transition zone are determined by linear interpolation of corresponding values for titanium carbide and (Ni-Cr) binder (that is proportional to the local concentration of titanium and carbon): where P is a considered parameter (Young modulus, Poisson ratio, yield stress, tensile or compressive strength and so on), C is a local concentration of corresponding constituent of the composite ( 1 ).The components of metal-ceramic composites have very different rheological properties (high-strength brittle particles of refractory compound and elastic-plastic metallic binder).To simulate the processes of deformation and fracture of such complex systems by MCA method the developed two-dimensional model of elastic-plastic interaction of cellular automata is used.An incremental theory of plasticity of isotropic medium with von Mises plasticity criterion and the plane stress approximation are used to model deformation of metallic binder on the mesoscopical scale.Elastic constants of the material and the diagram of uniaxial loading are used as input parameters for the model of interaction of cellular automata.These parameters determine mechanical response function of movable cellular automaton.To simulate the elastic-plastic metallic binder of the composite the parameters of mechanical response of movable cellular automaton conforming to the mechanical properties of nickel-chromium alloy were chosen.The response function of automaton modeling (Ni-Cr) is a "-" diagram with a linear hardening (curve 1 in Fig. 11).This diagram was obtained by approximation of the experimental diagrams for uniaxially compressed macroscopical samples of the alloy.Mechanical properties of automata simulating a high-strength brittle inclusion correspond to idealized properties of real particles of titanium carbide in ceramic phase (linear-elastic behavior up to failure, curve 2 in Fig. 11).An important role in building the mechanical model of composites belongs to determining the strength criterion (fracture criterion) for the system components (TiC and Ni-Cr) and the parameters of the chosen criterion.It is well-known that fracture is a fundamentally brittle and extremely localized phenomenon.Since the physical mechanisms of fracture are concerned with break of interatomic bonds and spatial separation of atomic layers, the criterion for material failure, in contrast to the criterion of plasticity, cannot be determined solely by the shear stresses and should take into account the effect of hydrostatic pressure.Therefore a two-parameter criterion of Drucker-Prager (36) was applied as one for interelement bond failure (fracture criterion) in the proposed model.It is necessary to note that for each of the composite components there are reliable published data about only one of the parameters of the criterion (36) 2. TiC: TiC c  =10000 MPa (tabulated value for ceramic phase), TiC t  =2000 MPa (the tensile strength was estimated under the assumption of the absence of significant defects in TiC particles), 5 As discussed above, the service characteristics of metal-ceramic composite (as well as of other composite materials) are largely provided by the adhesion of the binder to the reinforcing refractory inclusions.Depending on applied approach to the description of interphase boundaries the parameters of adhesive strength of dissimilar cellular automata (e.g., automata modeling TiC and Ni-Cr) were assigned in different ways: 1.When using the model of "narrow" interphase boundaries, the parameters bound  is equal to NiCr t  : Compressive strength for this pair is determined in a similar way:

Simulation results
A three-point bending test of the model samples of metal-ceramic composite TiC-(Ni-Cr) was simulated.Since the main research task was to analyze the changes in "macroscopical" response of the composite due to variation of the characteristics of interphase boundaries, the complex geometry of the mesoscopical TiC particles and features of particle size and spatial distributions were not taken into account.A composite with an idealized internal structure shown in Fig. 9 (the model of "narrow" interfaces) was considered.Note that TiC particles in real metal-ceramic composite may contain damages and microcracks.To eliminate the influence of this factor on the simulation results it was assumed that TiC particles in the model sample have a "monolithic" internal structure and do not contain significant defects.All interphase boundaries were assumed to be free of flaws as well.
The structure of the model set-up is shown in Fig. 12. Dynamic loading by 20 m cylindrical mandrel with constant velocity V load was applied to the sample of 24130 m in size.The values of loading velocity ranged from 0.2 m/s to 1 m/s.The value of the maximal resistance of the sample to bending and the corresponding displacement of the mandrel, the dynamics and pattern of the sample fracture as well as the fracture energy characterized by an area under the loading diagram were analyzed.

 
versus bending angle ".Here stress  is defined as a force of composite resistance to mandrel indentation referred to the area of the upper surface of the mandrel.The bending angle  is calculated as is shown in Fig. 13b.It is established that increasing the strength of "narrow" interphase boundaries leads to a twofold growth in strength of the composite and makes the value of the limit strain (strain corresponding to peak resistance) higher by an order of magnitude.Note that the falling (supercritical) part of the loading curves is associated with the initiation and development of main crack in the sample.), the value of composite strength reaches the saturation.The described stages are pronounced in the curve of interface strength depending on fracture energy (Fig. 14b).It can be seen that the fracture energy of composites with high-strength interphase boundaries increases by an order of magnitude in comparison with composites containing low-strength interfacial boundaries.


. Fracture energy was calculated as an area under the loading diagram expressed in the units "normalized stress -bending angle" (Fig. 13a).
The analysis of simulation results showed that staging of the change in integral parameters of the composite response with increase in strength of interphase boundaries ( bound t  ) is associated with a change in the nature of fracture.Fig. 15 shows the main stages of fracture in composites with different values of the strength of interphase boundaries.In the case of low values of interface strength bound t  (smaller than the strength of metallic binder NiCr t  ), from the early stages of loading (at small applied strain) TiC particles are detached (debonded) from the binder producing the extended splitted sites on the interphase boundaries, which are under tensile stresses (see top row in Increase in the applied strain joints these separate sites due to cracks passing through the binder.For high values of bound t  , under approaching of the strength of reinforcing particles the crack is nucleated, and it grows through the matrix bypassing the interphase boundaries, while interface zones remain virtually undamaged (see bottom row in Fig. 15).In the "intermediate" range of interface strength ( bound t  is equal to or a bit higher than the value of binder strength) the character of composite failure changes from an "interfacial" type to a "binder" one.It is clear that crack initiation in the plastic binder takes place under high applied strains, and crack growth is rather slow.As a consequence, fracture energy of composites with high-strength interfacial boundaries increases by an order of magnitude compared to composites with low-strength interfaces (Fig. 14b).Thus, the strength of interfaces between reinforcing inclusions and plastic binder is a critical factor, which determines a number of service characteristics in metal-ceramic composites, such as strength, critical strain, crack growth resistance and others.The strength of the interface is an integral parameter whose value is determined by a number of features of internal structure and geometry of the given boundary.One of the key factors determining the strength of metal-ceramic composites is the width of the transition zone (zone of variable chemical composition) H if , which determines gradients of concentration of chemical elements and of local mechanical properties.As shown in [52,55], shear stresses at the interfaces between elastic-brittle reinforcing particles and plastic matrices may exceed stresses in the matrix up to 5-10 times.This is due to a significant difference in the rheological characteristics of the interacting components of the composite.It is clear that high stress concentration determines a rapid debonding of particles from the binder in case of weak interphase boundaries (top row in Fig. 15).At the same time, wide transition zones at the interface can reduce the magnitude of the gradient of mechanical properties and, thus, decrease the stress gradient.That "smearing" of stress concentrator can lead to a substantial increase in the strength and deformation characteristics of the interfaces.This effect can be demonstrated by the given simulation results on three-point bending test of the metal-ceramic composite with "wide" interphase boundaries.A composite with the interphase boundary of 0.8 m wide was considered (the structure of the transition zone is shown in Fig. 10).Rheological and strength characteristics of movable cellular automata modelling the transition zone areas were determined on the basis of phase mixture model ( 43) using linear interpolation of corresponding values (parameters) for the titanium carbide and (Ni-Cr)-alloy.Fig. 16 shows the load-displacement diagrams for three-point bending test of the model sample of considered metal-ceramic composite (curve 1) and composites with "narrow" interphase boundaries (curves 2 and 3).It can be seen that a sufficiently wide transition zones provide high values of strength, ultimate strain and fracture energy of the material.The values of these parameters are close to the maximum corresponding to the metal-ceramic composites with high-strength "narrow" boundaries (Fig. 14).Note that for the composite with wide interfaces, fracture occurs by means of initiation and propagation of cracks in the binder.As a rule, the crack envelopes the transition zones in this case.The foregoing shows that the broad interphase boundaries with a gradual change in mechanical properties possess high values of the effective strength and have a significant influence on integral ("macroscopical") mechanical properties of metal-ceramic composite.As follows from analysis of stress distribution in model composites, high effective strength of wide interphase boundaries is a result of a significant expansion of the area of internal stress decreasing from high values in the reinforcing particles (which are stress concentrators in the composite [55]) to significantly lower values in the plastic binder (Fig. 17).This leads to two-three-fold decrease in the value of stress gradients for wide interphase boundaries as compared to "narrow" interfaces.As a consequence of a smoother distribution of internal stresses, the integral ("macroscopical") deformation and strength characteristics of metal-ceramic composite significantly increase, and fracture mechanism changes from damages and crack formation at the interfaces to the crack propagation in the plastic binder.Analysis of simulation results indicate that increase of the width of transition zone at the interphase boundaries in metalceramic composites can provide an increase in the integral parameters of mechanical response of materials (strength, ultimate strain and fracture energy).Fig. 18 shows dependencies of the sample strength (Fig. 18a) and fracture energy (Fig. 18b) in model metal-ceramic composite depending on width of transition zone.It is seen that presence of even narrow mesoscale transition zone (~0.3 m) between the TiC particles and (Ni-Cr)-alloy increases specimen strength and the fracture energy in 1.5 and 5 times respectively.A further widening of the interface leads to decrease of the rate of increase of integral parameters of the composite response.And when the value of H if becomes comparable to the average distance between TiC particles (H if >1 m) changing of strength and fracture energy of the composite becomes insignificant.

CONCLUSIONS
he group of particle-based numerical methods belonging to the concept of discrete elements (DE) is an extensively used and efficient tool to study complex processes of deformation and fracture of heterogeneous solids under complex loading conditions.The main advantage of these methods compared to conventional numerical methods in continuum mechanics (FEM, FDM and other) is a capability of direct modeling of numerous fracture accompanied by sliding and repacking of fragments, dilatancy of the medium and so on.Nevertheless, until recently, capabilities of this methods were limited by description of brittle and granular materials and media.A new approach which makes possible fundamental extension of the application field of DE-based methods to elasticplastic and visco-elastic-plastic solids is proposed in the paper.This approach is based on the idea about building associations between the components of stress/strain tensor of the local volume and the inter-element forces/displacements.The proposed associating allows one to rewrite relations of the applied model of plasticity (which are conventionally written in terms of stress/strain tensor components) in terms of forces and displacements or their increments.In particular, implementation of plastic flow theory with von Mises yield criterion within the discrete element concept is described.The movable cellular automaton (MCA) method [22][23][24] combining formalisms of cellular automaton methods and discrete element method was used as a numerical technique to implement this model of plasticity.Note that the developed approach provides potential possibility to realize various models of elastoplasticity or viscoelastoplasticity (including Dricker-Prager, Nikolaevskiy and other dilatant plasticity models) in the framework of "conventional" particle methods.Furthermore, it allows one to get isotropic (independent of the packing of elements) deformation pattern even on regularly packed particle ensembles, that is a fundamental problem in conventional formalism using pair-wise interaction potentials.Another important advantage of the developed formalism of discrete element interaction is a possibility to directly apply complex multiparametric fracture criteria (Drucker-Prager, Mohr-Coulomb, etc) as criteria of interelement bond breakage.The use of these criteria is very important for correct modeling of fracture of complex heterogeneous materials of various nature.
To demonstrate the capabilities of the approach some aspects of the problem of modeling of deformation and fracture of metal-ceramic composites are considered.On this example the prospects of application of particle-based numerical methods (in particular, of movable cellular automaton method) to study the influence of the features of internal structure of heterogeneous materials on their mechanical properties (strength, fracture toughness, etc) and fracture pattern are shown.At the present time described models of interaction of movable cellular automata belonging to the concept of discrete elements are approved and widely used to study fracture-related problems at different scales from nanoscopic to macroscopic one, whose investigation by conventional numerical methods of continuum mechanics is difficult.The problems of this type include physical and mechanical processes in contact patches of technical and natural frictional pairs [24,[56][57][58], multiple (quasi-viscous) fracture of porous ceramics or composite coatings [25,59,60], dynamics of damage accumulation and dilatancy in active fault zones [61,62] and so on.

T APPENDIX A
et's assume that under consideration of isotropic linearly elastic medium the expressions for specific forces of normal (central) and shear (tangential) response of discrete element i to mechanical action by neighbouring element j are defined on the basis of reformulation of Hooke's law for diagonal and off-diagonal components of stress tensor (19): Explanation of meaning of variables in (A1) can be found in Sections 2-4 (Main parameters of interaction of discrete elements -Discrete element interaction for modeling consolidated elastic-plastic medium).The proposed expressions (A1) are correct as they provide fulfillment of Hooke's law (19) for average stress and strain tensors in the volume of discrete element i.Actually, substitution of relations (A1) in the expression (16) for ==x leads to the following relationship: To understand the meaning of the first contribution in (A2) the notion of average strain tensor in the volume of discrete element i has to be introduced.Considering relative normal and shear displacements of interacting elements i and j as components of the vector of relative displacement the expression for average strains i   in the volume of element i can be formally written by analogy with i   .In considered two-dimensional problem statement it has the following form: According to (A3) the first contribution in the right part of (A2) has the meaning of corresponding average strain tensor component: To understand the meaning of the sum  ) has the form: where the identity was applied (the Kronecker delta   and the Einstein summation convention are employed here); r  is projection of a radius-vector to corresponding axis of laboratory system of coordinates.The volume integral in (A5) can be rewritten as a surface integral by applying the Gauss divergence theorem: where Si is the surface of element i; n  is the unit outward normal to infinitesimal surface area dS.The element i is loaded by forces acting at discrete interaction surfaces S ij (10).So, its total surface S i can be represented as a set of S ij , and surface integral can be replaced by the following sum: x n d S Sq where ij,  is defined as shown in the Ni is total number of interaction surfaces ( ).In accordance with (A7) the sum in the second contribution in the right part of (A2) is: So, the expression (A2) can be rewritten in the form: which is Hooke's law for average stress and strain tensor components i xx  and i xx  .
Corresponding expressions for other components of average stress tensor can be derived in analogous fashion: Relations (A9)-(A11) are derived on the basis of assumption that total surface S i of the element i is "occupied" by interacting neighbours j.However, if a part of surface S i of automaton/element i is free (i.e.neighbouring elements j "occupy" only part of Si), corresponding unoccupied surfaces Sik can be formally considered as surfaces of interaction with virtual neighbours having zero stiffness.In that case specific forces at unoccupied surfaces S ik has to be assigned as zero: ij=ij=0 (this doesn't hold true for corresponding pair strains: i(j)0 and i(j)0).Therefore equalities (A9)-(A11) are correct for discrete elements with partially free surface as well.In accordance with above mentioned the number of real neighbors of element i ( real i N ) can be used in expression ( 16)-(16) when calculating average stress tensor components.At the same time a total number of neighbours (including virtual ones) have to be used in expression (A3).Thus, average stress and strain tensor components are connected by Hooke's law when expressions (A1) are used for assigning central (normal) and tangential (shear) response of particle to mechanical impact of neighbours.This result confirms correctness of the proposed approach to description of interaction of discrete elements (or movable cellular automata) simulating locally isotropic linearly elastic medium.

APPENDIX B roposed expressions
for scaling specific central and tangential forces of response of the element/automaton i to the impact of the neighbor j provide a fulfillment of radial return algorithm of Wilkins ( 25)-( 26) for average stresses in the volume of the element i.
Actually, substitution of relations (B1) in the expression (16) for ==x leads to the following equation: According to (16) the sum 2 , , , 1 1 1 cos cos sin in the first contribution in the right part of (B2) has a meaning of corresponding average stress tensor component before stress returning procedure ("elastic" stress i xx  ).
In accordance with (A7)-(A8) the second sum   Thus, the proposed method (expressions (B1) and ( 26)) of scaling (return) of specific forces in pairs of discrete elements (or movable cellular automata) rigorously match radial return algorithm by Wilkins.This makes possible correct simulation of plastic deformation of consolidated solids under mechanical loading by the ensemble of discrete elements.

APPENDIX C
lgorithm of Wilkins for scaling local stresses in plane stress approximation (two-dimensional problem statement) is more complex than in general 3D problem or in plane strain state approximation.It is realized by iterative method [34].When using a single iteration only, the procedure of scaling of stress tensor components can be represented as the following sequence of operations done at each step of integration of motion equations.1.At the current time step (n+1) elastic problem is solved in incremental fashion.New values of stresses   and strains   are results of this solution.Calculated stresses can be further adjusted while strain values are "frozen".
2. New values of invariants of stress/strain tensor are calculated.Yield condition is checked.If the local stress state (point in the stress space) is inside the limit surface ( int pl    ), the procedure finishes.If the condition of crossing the yield surface is met, the value of stress scaling coefficient M is calculated.On of the features of stress scaling in plane stress state approximation is that the value of parameter M is a solution of the fourth-order polynomial equation [34]:   is mean stress value after stress scaling procedure (note that in contrast to 3D problem the procedure of stress correction in plane stress approximation is not accompanied by the persistence of mean stress because its value after the solution of elastic problem is overestimated [34]), index "n" indicates the value of corresponding variable referred to the end of the previous time step n.Implementation of this procedure within the framework of discrete element concept is concerned with scaling of element response forces  ij and  ij to satisfy the expressions (C3) for average stress tensor components in the volume of the element i.By the analogy with elastic problem the expressions for correction of ij and ij are obtained by direct reformulation of relationships for   and   in (C3): Note that in the case of multiple iterations of the procedure of stress/force scaling the relationships used have more cumbersome form [34], however their matter is the same.

Figure 1 :
Figure 1: Schematic presentation of switching between linked (at the left) and unlinked (at the right) states of the pair of discrete elements i and j.


of discrete elements i and j in the pair (in the general case  have opposite signs (in contrast to shear angles), as reflected in the expression(9).

Figure 3 :
Figure 3: An example of bending of the pair of linked elements i and j for the simple case:  i =  j , q ij = q ji and tan g

Figure 4 :
Figure 4: Examples of regular packing of elements in 2D problem statement: close packing (a) and square packing (b).
range for brittle materials) the value of S ij can rely constant ( 0 ij ij S S  ) ), sign «» implies use «+» in case of  = x and use «-» in case of  = y.Components i xz  and i yz  are identically zero.Definition of i zz  depends on constitutive equations of considered medium.Note that values of i xy  and i yx  coincide only in static equilibrium state of ensemble of discrete elements, while they can slightly differ at the stage of establishing static equilibrium.Therefore their mean value ( in the proposed model (hereinafter it is called as i xy  ).

Figure 5 :
Figure 5: An example of definition of angle  ij, between line connecting mass centers of discrete elements in the pair i-j and -axis of laboratory system of coordinates (=X is considered here).The center of coordinate system is translated to the mass center of the element i.

Figure 6 :
Figure 6: Schematic representation of functioning of radial return algorithm of Wilkins.Here  el is stress intensity after elastic problem solution at the current time step,  pl is a current value of yield stress (taking into account hardening),  y is an "initial" value of yield stress.


is stress intensity value, which results from solution of elastic problem at the current time step (n+1); is stress intensity value at the end of the previous time step n (after realization of stress return procedure if necessary).Note, that following the idea of Wilkins' algorithm the value of   int 1 i n

Figure 7 :
Figure 7: Instantaneous coordinate system XY associated with the current spatial orientation of the pair i-j.
ii. Calculation of current values of i xx  , i yy  , i xy  , i zz  and int i  with use of expressions (31).

used in ( 21 )
to get estimated values of  ij and corresponding estimated value of   1 predictor cycle.The number of predictor cycles is defined by assigned tolerance of mean stress change during one cycle (normally the solution converges very fast and 2 predictor cycles are enough).
of compressive strength ( c ) of the pair bond to tensile strength ( t ), int ij  and ij mean  are corresponding invariants of stress tensor ij      .Variables int ij  and ij mean  are calculated by analogy with (17)-(18).

minW
is the value of plastic work to make pair totally linked at threshold normal pressure  min ; max W  is the value of plastic work to make pair totally linked at normal pressure  max (work in(39) has to be considered in specific units (per unit of volume).Parameters min, max, min W  and max W  in criteria (

Figure 9 :
Structural model of metal-ceramic composite: (a) general view and (b) enlarged area inside of grey rectangle in (a).Main constituents of the structure are: nickel-chromium plastic binder (1) and high-strength brittle TiC particles (2).

Figure 10 :
Figure 10: The internal structure of a model metal-ceramic composite with a "wide" transition zone (~0.8 m in the present example) between TiC particles (shown in black) and the Ni-Cr binder (shown in light gray).Shades of gray color of varying intensity mark areas of the transition zone with appropriate content of TiC.
Tabulated data are the values of the compressive strength ( TiC c  ) for TiC and tensile strength ( NiCr t  ) for (Ni-Cr)-alloy.Values of the second parameters ( NiCr c  and TiC t  ) are unknown or estimated very roughly and uncertainly.Therefore, estimates of these missing parameters of fracture criterion were used in the model: 1. (Ni-Cr)-alloy: NiCr t  =700 MPa (tabulated value), t  and bound a are the only characteristics of the model interface.In the simulations carried out, the parameter bound t  was considered as a variable (assigned) value and varied within the limits between 0.5 NiCr t  (imitation of weak interfaces) and 2.5 NiCr t  (extremely strong interfaces).The parameter bound a in the criterion (36) was assumed to be 3 in all calculations.2. When using the model of "wide" interphase boundaries, the values of strength parameters t  and c  for the transition layer with certain local concentrations of components ( TiC C and NiCr C ) were determined on the basis of phase mixture model using the relation (1).For a pair of dissimilar cellular automata (characterized by different values of TiC C ) parameters c  and t  of bond fracture criterion (36) were defined as the minimal values corresponding to the strength characteristics of the areas modeled by these automata.For example, for a pair of dissimilar automata characterized by concentrations of TiC 1 TiC C  (automaton with pure TiC properties) and 0 TiC C  (automaton with pure (Ni-Cr)-alloy properties) the value of tensile strength for the pair t

Figure 12 :
Figure 12: Scheme of the simulation of three-point bending test of the model composite sample.The simulation results showed that the change in interface strength ( bound t  ) leads to a purposeful change in the integral characteristics of the composite response of 2-10 times.Fig. 13a shows diagrams of dynamic loading of simulated composite sample with different values of adhesive strength of the metallic binder to carbide particles (loading velocity in these calculations was V load = 0.4 m/s).The diagrams are plotted in the terms "normalized stress NiCr t

Figure 13 :
Figure 13: (a) The load-displacement diagrams for the model samples of metal-ceramic composite with different strength ( bound t  ) of

Figure 14 :
Strength (a) and fracture energy (b) in model metal-ceramic composite depending on interphase boundary strength bound t  .Values of interphase boundary strength and "macroscopical" composite strength are normalized by tensile strength of the binder NiCr t

Figure 15 :
Figure 15: Dynamics of fracture in metal-ceramic composite specimens characterized by different values of the strength of interface boundaries bound t  : 0.7 NiCr t 

Figure 16 :
Figure 16: The load-displacement diagrams for three-point bending test of model samples of metal-ceramic composite with "wide" (curve 1) and "narrow" (curves 2 and 3) interphase boundaries (V load = 0.4 m/s).Curves 2 and 3 correspond to composites with different values of strength ( bound t  ) of "narrow" interfaces: 1.0 NiCr t 

Figure 17 :
Figure 17: Stress intensity distribution in the central part of the samples of metal-ceramic composite with "narrow" (a) and "wide" (b) interphase boundaries.Three-point bending test.Stress distributions correspond to the bending angle  = 3.5

Figure 18 :
Strength (a) and fracture energy (b) in model metal-ceramic composite depending on width of transition zone.Value of composite strength is normalized by tensile strength of the binder NiCr t  .
the second contribution in (A2) the procedure of homogenization of unit second-rank tensor 1 considered.General expression for average value i e  of unit tensor component i e  in the volume of discrete element i ( part of (B2) is a diagonal component i xx e =1 of unit second-rank tensor.Therefore, the expression (B2) can be rewritten in the form: that use of expressions (B1) along with(26) provides for rigorous satisfaction of the necessary conditions

5 .
general, this equation has several different roots with different signs.The physical situation corresponds to the positive root that is less than 1 and closest to 1. General analytical solution of (C1) does not exist and therefore Wilkins proposed to use Newton's iteration method with int pl M    as an initial value for iterative procedure.3. Calculated coefficient M is then used to calculate irreversible (plastic) increment of Z-component of strain tensor ( variable szz in (C1)-(C2) corresponds to the solution of elastic problem at the current time step.4. Found parameters M and pl zz  are then used for scaling of szz and stress tensor components xx, yy, xy: marks scaled values of stresses.Note that original expressions in[34] bind scaling of diagonal stresses to correction (adjustment) of the values of plastic increments of strain tensor and volume change during the time step:   .Expressions (C3) are obtained by simple mathematical transformation of original relationships into the notation in terms of stresses: The total increment of Z-component of strain tensor (zz) is calculated: of strain increment, mean plastic) increment of Z-component of average strain tensor in the volume of the element i.It is easy to demonstrate that scaling of element/automaton response forces  ij and  ij in accordance with (C5) rigorously satisfies the expressions (C3) for i xx