Microsoft Word - numero_49_art_22_2485 N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 212 Focused on Russian mechanics contributions for Structural Integrity Algorithms for calculation damage processes Nikolay G. Burago Ishlinsky Institute for problems in mechanics of the RAS, Moscow, Russia burago@ipmnet.ru, http://orcid.org/0000-0002-1806-9386 Ilia S. Nikitin, Alexander D. Nikitin, Boris A. Stratula Institute of Computer Aided Design of RAS, Moscow, Russia i_nikitin@list.ru, http://orcid.org/0000-0003-3499-6910 nikitin_alex@bk.ru, http://orcid.org/0000-0002-2916-758X stratula@matway.net ABSTRACT. The paper reviews the existing approaches to calculating the destruction of solids. The main attention is paid to algorithms using a unified approach to the calculation of deformation both for nondestructive and for the destroyed states of the material. The thermodynamic derivation of constitutive relations for solids with elastic, viscous and plastic properties accounting possible destruction is presented. Explicit and implicit non-matrix algorithms for calculating the evolution of deformation and fracture development are presented. Implicit schemes are implemented using iterations of the conjugate gradient method, with the calculation of each iteration exactly coinciding with the calculation of the time step for two-layer explicit schemes. Therefore the solution algorithms are very simple. The results of solving typical problems of destruction of solid deformable bodies for slow (quasistatic) and fast (dynamic) deformation processes are presented. Recommendations are given for modeling the processes of destruction and ensuring the reliability of numerical solutions. KEYWORDS. Fracture; Damage; Elasticity; Plasticity; Thermodynamics; Softening; Grid methods; Through calculations. Citation: Burago, N.G., Nikitin, I.S., Nikitin, A.D., Stratula, B.A., Algorithms for calculation damage processes, Frattura ed Integrità Strutturale, 49 (2019) 212-224. Received: 15.04.2019 Accepted: 22.05.2019 Published: 01.07.2019 Copyright: © 2019 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. INTRODUCTION umerical computations for destruction processes have been conducted since computers and numerical algorithms appeared, i.e. more than 50 years. There are a lot of various approaches to obtain answers to such questions as which types of destruction processes are possible, when and where microdamages are going to originate, how long artificial and natural constructions can perform satisfying operational and safe conditions. In cases when a large margin of N http://www.gruppofrattura.it/VA/49/2485.mp4 N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 213 safety is required, for example in operation of buildings, it is enough to identify places with the highest possible stress levels and to take care that such stresses are several times smaller than the destructive ones for building materials. Such recommendations are formulated in the so-called theories of strength in the theory of resistance of materials. The next level of studying the processes of destruction is the study of the distribution of individual cracks to predict the possibility and time to continue the safe operation of structures in cases where the mentioned defects have already been found. Such problems are solved by analytical and numerical-analytical methods of brittle fracture mechanics [1]. Throughout the history of the application of numerical methods in this direction a very large amount of research has been performed. It is to admit that even for a single crack, the formulation of its propagation conditions and taking into account conditions (possibly contact) on newly formed free surfaces during its development is a very difficult mathematical task, which is practically impossible to implement in the case of multiple cracks. The number of subdomains between cracks is catastrophically increasing and their location is unpredictable in advance. It seems that the way to get out of this difficulty for fracture mechanics as well as for gas dynamics with multiple shock waves is to use through calculation methods with trapping narrow zones of large solution gradients. It may be done using zones for deformations localization that simulate cracks in solids, shock wave zones and contact discontinuities in gas dynamics. For the first time, through calculation methods for fracture modeling were described in Maenchen-Sack’s paper half a century ago [2]. However, for a long time, the insufficient performance of computers did not make it possible to effectively use such methods of through calculations. Now, when the performance of computers increased up by 3-4 orders of magnitude, the through calculation methods begin to dominate and allow not only to significantly expand the range of actual tasks to be solved but also to simplify the solution processes. Maenchen and Sack in their work considered the case in which material in an infinitely small volume was destroyed instantly when a certain criterion of destruction was fulfilled. The destruction was mathematically expressed in the replacement of the usual bond of stresses and strains in the framework of the theory of elasticity and plasticity by a bond describing the behavior of the destroyed material. In the principle axes of the stress tensor a new bond ensured the absence of resistance to tensile strength while maintaining the resistance to compression. It is important to notice that the solution algorithms, regardless of the loading rate beyond the yield point, are stepwise in time since the properties of materials depend on the loading history. In boundary value problems in the formulation of Maenchen and Sack, an important property of positive definiteness of the linearized operator of the problem is preserved for increments of the sought functions at each time step. Here, the criterion of correctness of Hadamard boundary value problems is satisfied [3]. This criterion ensures the existence, uniqueness and continuous dependence of the solution on the input data. Unlike the Maenchen-Sack approach (brittle destruction), in reality materials lose the ability to resist deformation not instantaneously but gradually so there is a long stage of softening. In the case of softening, step-by-step operators for solving the problems of elastoplasticity in increments use explicitly or implicitly the so-called tangential modules of elasticity (stress derivatives by strain). When softening, the tangent modules of elasticity become negative. Because of this at the time step the operators of the elastic-plasticity problems lose the property of positive definiteness, so the boundary problem becomes incorrect according to Hadamard, the numerical solutions become physically meaningless and the calculation ends abnormally. In continuum mechanics this phenomenon is known as a violation of the material resistance criterion according to the Drucker [4]. Thus classical theories of elastoplasticity are not suitable for describing gradual softening. A way out of this difficulty was found in [5, 6] where the experimentally observed drop in stresses with increasing deformations was due to degradation (decrease) of elastic modules and yield strength due to an increase in the density of microcracks, called damage. This meant that the areas of weakening did not reflect the properties of the material in alternating stresses - elastoplastic deformations, but were due to external causes (the growth of microdamages). So the use of negative tangent moduli of elasticity which violates the correctness of boundary value problems is not required. The introduction of damage allowed computing the processes of softening without violating the conditions of Hadamard and Drucker. To implement damage models in through calculation algorithms, the enhanced formulation of the constitutive relations for damage processes in deformable solids is required, that is suitable for describing the behavior of both the original intact and destroyed materials including the transition process of softening. Variants of the such constitutive relations may be found in [5, 7 – 14]. In the problems of calculating the fatigue failure of structural elements along with the proposing of the multi-axial fracture criteria under cyclic loading (see the review [15]) also proposed models that take into account the damage accumulation process with an increase in the number of loading cycles [16]. For example the corresponding differential equations for the damage function are formulated in [17]. The accuracy of through calculation methods for damage modeling may be significantly improved by the use of adaptive movable computational grids in order to minimize approximation errors in areas of large gradients of solutions (reviews and descriptions are in [18 – 20]). N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 214 In this paper we consider the above components of the formulation and numerical solution of problems on the destruction of solids, show how to implement the corresponding algorithms and give examples of solutions to typical problems on the destruction of solids. STATEMENT OF A GENERAL PROBLEM he system of equations describing the behavior of a thermoelastoplastic damageable medium is used here in the variant described in [9]. The system of equations contains the laws of conservation of mass, momentum and energy, as well as the kinematic relations: 1 0 1 0 0 0 1 1 ( ) ( ) 2 2 , T T T d d dU r dt dt dt d d dt dt                                            u e I e q F x ε I F F e L L ε x e ε L L ε L u u (1) and also the constitutive relations that deserve a more detailed consideration given below. The following traditional notation is used:  is density, u is a velocity, t is a time, x is an Euler radius vector (actual configuration), 0x is a Lagrangian radius vector (initial configuration), F is a deformation gradient, L is a velocity gradient, ε is an Almansi strain tensor, e is an Eulerian strain rate tensor, σ is a Cauchy stress tensor, U is an internal energy per mass unit, q is a heat flux vector, T is a temperature, r is a mass heat source, /d dt is a material time derivative,  is a spatial differentiation operator in the current configuration and I is a unit tensor. Constitutive equations represent the relations between the characteristics of the state of an infinitely small volume of material, imposed by the laws of thermodynamics. Let’s arrange a minimal set of mutually independent state parameters of an infinitely small volume of the continuum: 0 0 0 0 ddTT T dt dt       χ ε χ e , where 0( )p  χ ε are structural parameters namely 0pε is a plastic strain tensor and  is a damage and is defined further. They change the internal structure of the continuum, so they are responsible for the development of dislocations and microcracks. Marked with zeros are material tensors associated with spatial tensors by the ratios 0 0 0 1 T T T           ε F ε F e F e F σ F σ F The law of entropy increasing  from the first law of thermodynamics, which asserts the law of conservation of energy, and the second law of thermodynamics 0 d r dt T T           q whence the inequality of dissipation rate is 0 0 0 ,0 0 0t dT T D T dt T                                 σ e χ q ε χ Here U T   is free energy per mass unit. Free energy and dissipation rate are T N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 215   2 0 0 0 ( )) ( ) ( ) 2 p p TK ln T T                    ε ε ε ε 2 ( ) ( ) q p p p p kd D H k H k T T dt T                   e e where H is the Heaviside's function. It is assumed that the elastic components of the strain deviator are small compared to unity. The effect of thermal expansion is taken into account by a member with a coefficient  . The component of the dissipation rate, which is responsible for the plastic flow, is assumed to be a homogeneous function of the first order of the rate of plastic deformation, which corresponds to the case of an elastoplastic medium. Plastic deformation increases when active loading condition 0 0( ) 0p pT      ε ε e is satisfied. It is also assumed that the material is plastically incompressible i.e. the dissipation rate depends only on the plastic strain rate deviator, which is usually well performed for metals. The resistence of the medium is represented by the modules of elasticity and yield strength. In addition to temperature, deformation and plastic deformation it also depends on the additional structural parameter of the state  . This parameter is called damage. Following relations utilize the damage parameter: 0 ( )g   , 0 ( )KK K g  , 0 ( )p p pk k g  , where 0 and 0K are elastic modules and 0pk is a yield strength for intact material. The functions g , Kg and pg are decreasing from 1 to 0 while   . They provide a decrease in the resistance of the medium with an increase in the damage that occurs when the condition of destruction 0 0( ) 0pT      ε ε e is fulfilled. The kinetics of the destruction process is determined by the dependence of the dissipation rate on the growth rate of damage. Non-negative functions 0 , 0K , 0pk , k and qk depend on 0T ε and 0pε . So the constitutive relations take the following form: 0 0 0 1 2 ( ) ( ) ( ) ( ) q p p p p p p U T k T p T T p K ln T T d H k H k dt                                                               q σ I σ σ ε ε e e e (2) Boundary conditions are 0 0 ( ) 0 0 ( ) 0 0 ( ) un n un n u u T T n S t u S \ S t p S t u S \ S t p S t T T S \ S t q                                                    x u n x σ n n x u x σ n x x q n (3) where 1, 2  . Initial conditions are 0 0 00 0 0px V t T T             x x u u ε (4) Thus it is required to solve the initial-boundary problem for the system of equations (1), (2) with the boundary (3) and initial (4) conditions. N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 216 Appendix to section 2 The purpose of this appendix is to clarify the meaning of the introduction of the damage parameter. Let us consider the system of equations of the Prandtl-Reuss classical nonlinear elastic-plastic model: / : ( )pd dt     u σ σ E ε ε 0( ( )) /s d dt     ε x x x u / ( ) : ( ) 0p p p pd dt H     pε λ σ ε ε where ( ) ( ) 2Ts   A A A is a symmetrization operator. For increments of stresses and strains at a time step the following relations take place :t  σ E ε where fourth order tensor tE depends on full and plastic deformations as well as on the loading mode (active or unloading). Let’s rewrite the stress increments in the following form : ( )t s    σ E x For velocities at each time step we have the following equation / ( : ( ) )sd dt t      tu σ E u The correctness of this equation is determined either by the Hadamard condition 0tE   σ or by its physical equivalent known as Drucker condition 0t t  σ ε Destruction is accompanied by weakening 0tE  and leads in the framework of the considered classical theory to the loss of correctness of the problem (unlimited exponential increase in speed, loss of continuous dependence of the solution on the input data, instability). Let’s consider the phenomenon of destruction using the following one-dimensional problem of tensile rod. Rod’s material is elastic with Young’s modulus E . Let’s assume that there is a small area of weakened resistance (“fracture zone”) with a small Young's modulus 1E E in the middle of the rod. It may be shown by an elementary solution that there is a surge of deformations in the “fracture zone” (Fig. 1a) and the displacements (and velocities) change almost by jump (Fig. 1b). The question arises about possibility to construct a physically and mathematically correct model of an elastoplastic material that describes the appearance of zones of reduced material resistance similar to the “fracture zone” from the considered model elastic problem. The affirmative answer is obtained from the theories of damage [5, 6]. Indeed, with the introduction of the additional sought- for damage parameter  the system of elastoplasticity equations takes the following form: / ( ) : ( )pd dt E     u σ σ ε ε N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 217 0( ( )) /s d dt     ε x x x u / ( ) : ( ) 0p p p p pd dt H       pε λ σ ε ε / ( ) ( ) 0 ( )p pd dt H                 ε ε ε ε where H is the Heaviside's function. For increments of stress at the time step we get: : ( ) ( ) : ( ) 0p p            σ E ε ε E ε ε E E a) b) Figure 1: An exact solution for a stretchable rod with a weakened central part: qualitative graphs of the distribution of deformation (a) and displacement (b) along the rod. The drop in stresses regardless of the change in strain is provided by the second term in the above expression for stress increments. The coefficient 0E  is negative, which reflects the fact that the modules of elasticity ( ) 0E   and stress decrease with increasing damage. Fig. 2 may serve as illustration of this process. The positivity of the elastic modulus ensures the fulfillment of the necessary conditions for the correctness of the initial-boundary problems according to Hadamard and Drucker. It should be noted that, together with the elastic moduli, the yield strength also decreases with increasing damage. The left-hand image in Fig. 2 shows the experimentally observed dependence in stress-strain variables. The right-hand image in Fig. 2 shows the same process in the extended stress-strain-damage coordinates. One may see that for any fixed value of damage, there are no areas of weakening in the stress-strain diagrams. Figure 2: The effect of damage on the deformation diagrams. In the theory of damage the responsibility for softening is no longer on the condition of plasticity but is on the equation of damage. Мore of that the damage may take place independently of plastic deformation The softening or, more generally, the destruction is considered as the loss of the ability of the material to resist elastic deformation, expressed in a decrease in the values of elastic modules and yield strength due to breaking of elastic bonds or, equivalently, due to an increase in the N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 218 number of microcracks per unit volume. Destruction is described as a process independent of deformation (exactly independent). Therefore, the fracture characteristic called the damage parameter is related to deformation, temperature, and other state parameters only to the extent that they all participate in the general system of equations and initial-boundary conditions of thermomechanics. In fact, the destruction can happen without any deformation at all, for example due to chemical reactions or laser radiation and other external influences of a non-mechanical nature. It is important to note that the areas of weakening of experimental diagrams of material deformation cannot be used to determine the properties of elastoplastic materials since in these areas the deformation processes are controlled by increasing damage through the change of elasticity coefficients of the material and yield strength. The formulation of theories and reasoning were deliberately simplified in this commentary. In particular, temperature, hardening parameters and hardening dependence on strain rate were not considered; members related to large deformations were neglected. The goal was to more clearly describe the basic idea of damage theories. It has to be noticed that these simplifications are not used further as the basis for our numerical calculations. For describing the elastic properties of material the simplified form of Hooke's law is used, obtained under the assumption of the initial isotropy of the properties of the material: ( ) : 2 ( )p p         I I where the stiffness constants of Lame and the yield point, which determines the boundaries of the elastic behavior of material, depend on damage as follows: 1000 0e    , 10000e    , 10000p p e    . Values for the intact material are marked by index "0". We assume that the damage grows if the maximum principal deformation reaches the positive critical value. So the maximum principal deformation should be a tensile deformation. Destruction condition takes the form:  1/22 21 ( ) ( ) 4 0 2 d x y x y xy dF M                where the dimensionless scale factor is introduced max min max min min( , ) max[( ),( )] x yh h M x x y y    It represents the well known feature, inherent in the concentration of deformations in the tip of crack in the elastic material, which makes it possible to treat this local criterion as the approximation of the usual criterion of destroying the mechanics of brittle failure using of the coefficients of stress-strain concentration. The criterion of destruction with the scale factor ensures the convergence of the numerical values of the integral critical loads of destruction under the mesh refinement. Otherwise since the concentration of deformations is described increasingly better on the finer grid, the critical deformation of destruction is reached increasingly earlier, with the smaller values of the applied load, that tends to zero with grid refinement. Certainly, calculated critical loads depend on grid permission, but this dependence must demonstrate the convergence of calculated critical loads to a certain nonzero limiting value. Such behavior of nymerical model ensures the introduced above scale factor. The introduction of such coefficient is not the final solution of the problem of the convergence of the calculated critical loads of destruction, but it is possible that this is step to the side of the possible solution of described problem. It can be, someone will devise the best regularization of the criterion of destruction. Under the conditions of the absence of the information about the kinetics of damage, we accept the kinetic equation for the damage in the simplest form: 1000 ( )t dH F  where H is the Heaviside's function. Large magnitude coefficient in the right side is introduced for guaranteeing the rapid growth of damage in order to provide the lost of material resistance for several temporal steps. This makes it possible to N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 219 consider mathematical model of damage used here as regularized and, simultaneously, the simplified version of Maenchen- Sack model [2], in whom the stress-strained state due to the destruction is corrected instantly by jump. The instantaneous correction of stresses leads to the appearance in the numerical solution of the nonphysical oscillations and requires smoothing of sought functions is required. NUMERICAL METHOD he solution algorithm is based on the modification of an implicit finite element scheme built in [20] and implemented as part of the ASTRA application package. The main features of the algorithm are as follows. The initial equations of the problem including the constitutive and kinematic differential relations are applied in the integral variational form of Bubnov – Galerkin. A simple piecewise linear finite element approximation is introduced by spatial variables on a moving grid containing triangular and quadrilateral cells in the two-dimensional geometry and tetrahedra, prisms and parallelepipeds in three-dimensional geometry. Piecewise linear approximation is applied to all desired functions, including displacements, velocities, temperature, heat fluxes, deformations, plastic deformations, stresses, hardening parameters and damage, the discrete values of which are represented by nodal values. The set of basic functions contains displacement, velocity, temperature, plastic deformation and damage parameter. The rest of unknowns can be determined by basic functions using spacial differentiation and non-differential relations. Numerical integration points are located at the nodes of the grid, so the mass matrix is diagonal. At each time step the nonlinear terms of the equations are linearized by the Newton method with respect to small increments of the basic functions. To solve auxiliary linearized problems iterative method of conjugate gradients is used. It is implemented without matrix operations and working on each iteration as an explicit two-layer central-difference scheme. As the time step decreases, the number of iterations of the conjugate gradient method for solving linearized equations decreases (due to better initial approximation) and for time steps within the Courant constraint the implicit scheme works asymptotically as fast as usual explicit schemes. Under the least successful initial approximation no more than N iterations are required to determine the solution, where N is the number of unknowns. Since physical splitting is applied, the velocities, temperature, plastic deformation, damage and coordinates of the moving grid nodes are determined in turn and the number of discrete variables is not too large. The time step is chosen in implicit schemes from the accuracy conditions in order either to limit the maximum deformation increment to about one tenth of the deformation corresponding to the yield strength or, for elastic materials, to a value much less than one. In explicit schemes the usual stability conditions are applied for the cases of hyperbolic and parabolic equations. For preconditioning of the algebraic problems (multiplication of the system of equations by an approximate inverse matrix) diagonal approximations of stiffness matrices are used that corresponds to scaling of unknowns. This is quite enough to ensure the stability of conjugate gradients iterations, even if penalty terms are used to account contact or incompressibility conditions. The conditions of approximation, stability, and convergence of numerical solutions are justified on the level of simplified model equations. Therefore, even a stable numerical solution may contain unwanted non-physical oscillations. The wavelength and period of such oscillations coincides with the length of the steps of the computational grid in space and time. A large number of various approaches are known to eliminate such disturbances. Usually the introduction of explicit or scheme viscosities makes the solution monotonous. For problems of destruction such smoothing methods should be introduced very carefully, because the physical phenomenon of destruction is very sensitive to small disturbances. Too coarse smoothing procedure (such as averaging in the Lax method) or, conversely, its absence can strongly distort the solution. Localization of deformations does not take place if the non-physical viscosity of the scheme is large, or, conversely, the location of the fracture zones will be strongly and obviously dependent on the grid steps in time and space if the smoothing is insufficient. In the present work, the following method of monotonizing numerical solutions is recommended and utilized. At the end of each time step the function f that is to be monotonized is twice differentiated in each of the spatial independent variables. Smoothing process is performed in a coordinate-wise manner although the grid may contain cells with arbitrarily spaced edges. When using a piecewise linear approximation, the values of the second derivatives xxf are calculated as generalized solution from the variational equation: ( ) 0xxxx xx V ff f f dV x x        T N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 220 Since it is assumed that at the boundary points the first or second derivatives are negligible, in the equation for the second derivatives the integrals over the boundary of the solution domain are rejected. Then for each node of the grid the averaged value f of the monotonized function is calculated. This value also depends on the coordinate direction (to put it simply, it is the average of the function values at points of intersection of the border of the neighborhood of a node with a given coordinate line). Finally the old values f are substituted with new ones ( ) 2f f  only in those few grid nodes (points of non-monotonicity) in which the second derivative xxf changes the sign at the nodes of one edge. We want to emphasize an important fact that the old value is not replaced by the average, as in the Lax method, but is shifted only at half towards the average value. Unlike most other methods the described monotonization procedure does not work in all grid nodes, but changes the solution only at non-monotonic points. Acting more selectively it can even turn an unstable scheme into a stable and converging to the right solution one. It well regularizes standard two-layer central-difference explicit schemes not only for problems of the dynamics of solid deformable bodies, but also for problems of Euler mechanics of fluid and gas. Of course, when constructing such schemes, the elementary requirements for the description of dissipative processes at jumps must be taken into account i.e. destabilizing terms with negative viscosity coefficients found in the first differential approximations of such difference schemes must be balanced by adequate artificial viscosity. In this case this method of monotonization works fine as a supplementary one. Here the monotonization procedure is used in conjunction with an implicit scheme as a control of the monotony of the solution. The procedure does nothing if the solution is monotonous but immediately eliminates any non-physical oscillations if they occur. It should be mentioned that described monotonization procedure is not conservative. That is why additional conservation control and correction procedure is desirable. RESULTS elow we examine the task about the development of the zones of destruction in the stretchable flat square plates with the round and elliptical macro-pores or rigid inclusions. The two-dimensional spatial domain can be seen in Fig. 3-5. We use the grid of triangular finite elements is usual almost uniform grid, generated automatically. On the vertical boundaries of the square region of the solution horizontal displacements is relied by equal to zero. Lower boundary has zero vertical displacement. Upper boundary slowly moves up (with the speed that is much less than speed of sound in the material) and provides quasi-static extension of the body in question. Shearing stresses on the boundaries are relied by equal to zero. On the surface of the macro-pore external forces are absent. Rigid inclusions are filled with the absolutely rigid material, completely coupled with the material of the body in question. Finite elements in region of rigid inclusions have huge values of the elastic moduli and yield point (100 times more than for the plate material). Calculated diagrams of deformation depict the dependence of the stress averaged over the upper horizontal section on the monotonically increasing assigned vertical displacement of upper boundary. In the dimensionless coordinates in the solution domain are: 6.0 6.0x   и 6.0 6.0y   , the diameter of circular pore and rigid inclusion is equal to 1.0, the semiaxes ratio of elliptical pores and inclusion is equal to 0.5, angle of rotation is equal 30 . Young's modulus is equal to 1000, Poisson ratio is equal to 0.2, the speed of sound is equal to 1.0, the speed of the motion of upper boundary yv grows from zero for 0t  to the value 0.001 for 510t  . In the calculations according to the model of ideally-plastic material the yield point is equal to one. The deformation of destruction is equal to 0.02. The results for the elastic plate with the circular rigid inclusion are shown in figure 3. Black narrow zones correspond to the developing macrofissures, darkening answers the amount of maximum principal deformation. Right graph depicts the calculated diagram of deformation. The case with one circular pore is demonstrated in figure 4. A difference in the pictures of destruction for inclusion and pore is explained by the different nature of the concentration of deformations near the pore and near the rigid inclusion (Fig. 5). In the case of the pore the destruction begins from the points in the horizontal direction outermost from the center of the pore, and for the rigid inclusion destruction earlier begins at the points in the vertical direction outermost from the center of inclusion. Solution for the elliptical pore, oriented at the angle, in the elastic-plastic material it is shown in figure 6, where are depicted the zone of destruction (a), the calculated diagram of deformation (b), the distribution of plastic work and vertical B N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 221 displacement. The splash of the values of plastic work in the zone of destruction and an jump like change of the displacement is distinctly visible. a) b) Figure 3: Circular rigid inclusion in the elastic material: the zone of destruction and the distribution of the maximum principal deformations of (a) and the calculated diagram of deformation (b). a) b) Figure 4: Destruction of elastic material with the presence of the circular pore: the zone of destruction and the distribution of the stress y (a); calculated stress-strain diagram (b). a) b) Figure 5: Maximum principal deformations near the pore (a) and the rigid inclusion (b) in the elastic material N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 222 a) b) c) d) Figure 6: Destruction of ideal elastic-plastic material with the presence of one large oval pore: the zone of destruction (a); the calculated diagram of deformation (b); the distribution of plastic work (c); vertical displacement (d) In figures 7 and 8 the zones of destruction and the calculated diagrams of deformation for the ideal- elastic-plastic models with the large and small elliptical pores are compared. Result is obvious: with an equal quantity of pores the sample with small pores has greater strength. a) b) Figure 7: Destruction of ideal-elastic- plastic material with the group of the large oval pores, oriented at the angle 30o , in the absence the horizontal displacements of the lateral boundaries: the zone of destruction (a); the calculated diagram of deformation (b). In the parametric analyses the influence of the characteristics of discreteness is investigated. Depending on the size of time- spatial steps and on the form of mesh the local picture of destruction is different in details. However, the integral averaged characteristics of process (represented by diagrams of the deformation) are practically the same. N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 223 Below we present results for cracking of massif near drill-hole. The continuum model of stochastically inhomogeneous rock massif destruction was implemented for modeling the formation and development of fractal fracture structures in the vicinity of a cylindrical hole under the excessive internal hydrostatic pressure (Fig. 9). a) b) Figure 8: Destruction of ideal-elastic-plastic material with the group of the small oval pores, oriented at the angle 30o , in the absence the horizontal displacements of the lateral boundaries: the zone of destruction (a); the calculated diagram of deformation (b). a) b) Figure 9: Formation and development of fractal fracture structures in the vicinity of a cylindrical drill-hole under the excessive internal hydrostatic pressure. The range of random strength limit values 5% (a), and 20% (b) The criterion of fracture based on the limit values of principal deformation is used. The stochastic heterogeneity of the strength properties of the massif was provided by an artificially limited spread of the values of the parameters of the damage criterion using a random number generator. Stochasticity ensures the formation and development of multiple zones of fracture, forming fractal structures. The introduction of the initial spread of the strength properties of the medium is quite capable of reproducing the cracking of rocks with the formation of fractal fracture structures. With the same loading programs, the damage patterns are quite different due to difference in the level of stochasticity of the strength properties (Fig. 9-a, b). CONCLUSIONS umerical experiments have shown that fracture criteria for ultimate deformations for elastoplastic materials work better than criteria for ultimate stresses. The following properties of the theoretical model also play an important role to obtain localization zones in the form of narrow crack-like bands:  An abrupt decrease in material’s resistance with increasing damage;  The accuracy of calculation controlled by increments of deformations;  Accounting for inertia forces (natural regularization of the problem); N N. Burago et alii, Frattura ed Integrità Strutturale, 49 (2019) 212-224; DOI: 10.3221/IGF-ESIS.49.22 224  Minimum smoothing (especially in areas of destruction);  Resistance of the fractured material to compression (to ensure the nondegeneracy of the law of motion). Due to the above properties, the usual model of a damaged elastoplastic medium describes well the sharp localization of deformations along thin bands of large gradients of displacements and deformation and supports the convergence of the solution in the process of developing narrow fracture zones. The issues of bringing the model of a damaged elastoplastic medium in accordance with the data of physical experiments remain open and relevant. By using ASTRA application package as it was done here one is able to confidently carry out such calculations. The number of calculated examples, as well as the degree of detail of the description of the results can be easily increased. ACKNOWLEDGMENTS he work was supported by RSF, grant № 19-19-00705. REFERENCES [1] Cherepanov, G.P. (1979). Mechanics of Brittle Fracture, New York, McGraw Hill. [2] Maenchen, G., Sack, S. (1964). The TENSOR code, In: Methods in Computational Physics, v.3, Fundamental methods in Hydrodynamics, New York, Academic Press. [3] Hadamard, J. (1902). Sur les problèmes aux dérivées partielles et leur signification physique, Princeton University Bulletin, pp. 49-52. [4] Drucker, D.C. (1964). On the postulate of stability of material in the mechanics of continua, Mekhanika. Period. Sbornik Perevodov Inostr. Statei, 3, pp. 115-128. [5] Kachanov, L.M. (1986). Introduction to Continuum Damage Mechanics, Kluwer, Acad. Publ. [6] Rabotnov, Yu. N. (1959). Mekhanizm dlitelnogo razrusheniya, Sbornik “Voprosy prochnosti materialov I konstruktsii”, pp. 5-7. [7] Lemaitre, J. (1996). A Course on Damage Mechanics, 2nd ed., Berlin, Heidelberg, New York, Springer. [8] Fomin, V.M., Gulidov, A.I., Sapozhnikov, G.A. et al. (1999). Vysokoskorostnoe vzaimodeistvie tel, Novosibirsk, Izdatelstvo SO RAN. [9] Burago, N.G., Glushko A.I., Kovshov A.N. (2000). Termodinamicheskii metod polucheniya opredelyayuschih uravnenii dlya modelei sploshnykh sred, Izvestiya RAN. Mekhanika tverdogo tela, 6, pp. 4-15. [10] Bazant, Z.P. (2002). Reminiscences on four decades of struggle and progress in softening damage and size effect, Concr. J. (Japan Concr. Inst.), 40, pp. 16-28. [11] Kondauro,v V.I., Fortov, V.E. (2002). Osnovy termodinamiki kondensirovannoi sredy, Moscow, MFTI. [12] Krajcinovic, D. (1996). Damage Mechanics, Amsterdam, Elsevier Science. [13] Oliver, J. (1999). Continuum modeling of strong discontinuities in solid mechanics using damage models, Comput. Mech., 17, pp. 49-61. [14] Voyiadjis, G.Z., Kattan, P.I. (1999) Advances in Damage Mechanics: Metals and Metal Matrix Composites, Amsterdam, Elsevier. [15] Burago, N.G., Zhuravlev, A.B., Nikitin I.S. (2011). Models of multiaxial fatigue fracture and service life estimation of structural elements, Mechanics of Solids, 46, pp. 828-838. [16] Lemaitre, J., Desmorat, R. (2005). Engineering Damage Mechanics, Berlin, Heidelberg, Springer. [17] Marmi, A.K., Habraken, A.M., Duchene, L. (2009). Multiaxial fatigue damage modeling at macro scale of Ti6Al4V alloy, Int. J. of fatigue, 31, pp. 2031-2040. [18] Ortiz, M., Quigly, J.J. (1991). Adaptive mesh refinement in strain localization problems, Comput. Meth. Appl. Mech. Engng., 90, pp. 781-804. [19] Liseikin, V.D. (2010). Grid generation methods, Springer. [20] Burago, N.G., Nikitin, I.S., Yakushev, V.L. (2016). Hybrid numerical method for unsteady problems of continuum mechanics using arbitrary moving adaptive overlap grids, Computational Mathematics and Mathematical Physics, 56, pp. 1064-1074. T << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles true /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Error /CompatibilityLevel 1.4 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.0000 /ColorConversionStrategy /CMYK /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams false /MaxSubsetPct 100 /Optimize true /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments true /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile () /AlwaysEmbed [ true ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 300 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.50000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 300 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.50000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.50000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile () /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName () /PDFXTrapped /False /CreateJDFFile false /Description << /ARA /BGR /CHS /CHT /CZE /DAN /DEU /ESP /ETI /FRA /GRE /HEB /HRV (Za stvaranje Adobe PDF dokumenata najpogodnijih za visokokvalitetni ispis prije tiskanja koristite ove postavke. Stvoreni PDF dokumenti mogu se otvoriti Acrobat i Adobe Reader 5.0 i kasnijim verzijama.) /HUN /ITA /JPN /KOR /LTH /LVI /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken die zijn geoptimaliseerd voor prepress-afdrukken van hoge kwaliteit. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.) /NOR /POL /PTB /RUM /RUS /SKY /SLV /SUO /SVE /TUR /UKR /ENU (Use these settings to create Adobe PDF documents best suited for high-quality prepress printing. Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.) >> /Namespace [ (Adobe) (Common) (1.0) ] /OtherNamespaces [ << /AsReaderSpreads false /CropImagesToFrames true /ErrorControl /WarnAndContinue /FlattenerIgnoreSpreadOverrides false /IncludeGuidesGrids false /IncludeNonPrinting false /IncludeSlug false /Namespace [ (Adobe) (InDesign) (4.0) ] /OmitPlacedBitmaps false /OmitPlacedEPS false /OmitPlacedPDF false /SimulateOverprint /Legacy >> << /AddBleedMarks false /AddColorBars false /AddCropMarks false /AddPageInfo false /AddRegMarks false /ConvertColors /ConvertToCMYK /DestinationProfileName () /DestinationProfileSelector /DocumentCMYK /Downsample16BitImages true /FlattenerPreset << /PresetSelector /MediumResolution >> /FormElements false /GenerateStructure false /IncludeBookmarks false /IncludeHyperlinks false /IncludeInteractive false /IncludeLayers false /IncludeProfiles false /MultimediaHandling /UseObjectSettings /Namespace [ (Adobe) (CreativeSuite) (2.0) ] /PDFXOutputIntentProfileSelector /DocumentCMYK /PreserveEditing true /UntaggedCMYKHandling /LeaveUntagged /UntaggedRGBHandling /UseDocumentProfile /UseDocumentBleed false >> ] >> setdistillerparams << /HWResolution [2400 2400] /PageSize [612.000 792.000] >> setpagedevice