MODIFICATION OF IMPLICIT ALGORITHM FOR SOLVING A PROBLEM ON THE ELASTIC PLASTICITY OF BULK MATERIALS

A mathematical statement is given of the elastic-plastic behavior of isotropic bulk material using a classic model of Drucker-Prager. We have improved the return-mapping algorithm for solving numerically a problem on the mechanical state of bulk material. In order to solve a system of nonlinear equations by the Newton method, it is proposed, at each step of iterations, instead of finding the inverse matrix, to solve a system of algebraic equations, linearized by Newton, by applying the Gauss exclusion method. This makes it possible to reduce the number of arithmetic operations by about 3n2 (n is the dimensionality of SLAE) at each iteration step for each plastic finite element. We have tested the programming code developed on the high-level programming language Fortran on the example of a model material, characterized by the associative law of current, at different values of the angle of natural repose. Comparison of the obtained results of numerical experiment with the data received by applying the proprietary software revealed a deviation within 0.25–5.3 % depending on the desired magnitude


Introduction
Bulk carbon-graphite materials are widely used as high-temperature thermal insulators and heating elements in resistance furnaces of the electrode-manufacturing industry [1,2].To account for the dependence of physical properties of bulk materials on pressure in the course of numerical analysis of the thermal-electric state of high-temperature furnace equipment of electrode production [3,4], it is necessary to know pressure distribution in the layers of these materials.Failure to take this dependence into account may lead to significant errors in the results of numerical analysis during development of new equipment and when working out the norms of its operation.It is known that the information on pressure distribution in bulk materials in the form of the mean hydrostatic pressure can be obtained from the solution to a nonlinear problem on plasticity by the Drucker-Prager yield criterion [5,6].That is why the task to improve algorithms for solving the problem on mechanical behavior of bulk materials is a relevant issue.

Literature review and problem statement
An analysis of data from the scientific literature revealed that the Drucker-Prager model had become wide-spread in the numerical studies into elastic-plastic state of loose and brittle materials: -examining the strength of concrete with different composition [7][8][9]; -the process of rock destruction [10,11]; -a mechanism of compaction of pharmaceutical powder materials [12,13]; -the process of formation of metal-powder-like products [14][15][16][17].
Paper [7] describes the use of the extended Drucker-Prager yield model during numerical simulation of the behavior of reinforced concrete.Numerical experiments were performed employing the proprietary software ABAQUS.
Article [8] describes experimental studies into determining the limits of strength at compression of cylindrical concrete samples during uniaxial compression.The method is given for graphical determination of the Drucker-Prager criteria for concrete during its uniaxial loading.
In order to obtain the yield surface by Drucker-Prager for samples made of different grades of concrete, paper [9] experimentally determined coefficients of cohesion and the angles of inner friction.
Article [10] presented a model of mechanical behavior of rocks with different porosity, which is based on the combination of the Drucker-Prager boundary surface
Paper [11] considered a model of deformation of rocks that is a generalization of the Hill model of anisotropic plasticity, on the one hand, and the Drucker-Prager model, on the other hand.Underlying the model is the non-associative law of plastic current with strengthening for an anisotropic body, which takes into account the impact of volumetric stresses.The authors considered a rather general case of the combination of isotropic and translational strengthening.
Article [12] reports results of numerical research into mechanism of compaction of powder material into tablets.By using the commercial software ABAQUS, the authors performed numerical study in an axisymmetric statement.Solving an elastic-plastic problem employing the Drucker-Prager model made it possible to determine pressure distribution in the material in the form of the mean hydrostatic pressure, which is necessary for determining physical properties of the resulting material.
Paper [13] describes numerical and experimental studies into the pressing process of pharmaceutical powder materials.Numerical studies were performed using the Drucker-Prager model.The ABAQUS software was employed.
In [14], numerical studies are reported for the pressing process of aluminum oxide powder applying a modified Drucker-Prager model.Numerical experiment was performed using the ABAQUS software application.
The authors of [15], by employing the method of finite elements, carried out research into processes of plastic deformation of metallic powder materials.It is shown that the ABAQUS software package enables obtaining more accurate results for the Drucker-Prager model than when using the software DEFORM and ANSYS/LS-DYNA.
Papers [16,17] report physical and numerical experiments of the pressing process of a mixture of metallic powders.Numerical study is represented in an axisymmetric 3D statement using the method of finite elements.The calculations were carried out applying a modified Drucker-Prager model implemented in the software ABAQUS.
Certain shortcomings of the considered studies that address application of the Drucker-Prager model for loose and brittle materials include the following: -a lack of complete mathematical formulation of the problem and the solving algorithm; -there are formulae only to determine the yield criterion, equivalent deformations and stresses, results and analysis of numerical studies, etc.; -the calculations are carried out using the commercial software ABAQUS, employing which requires purchasing an appropriate license.
Thus, the studies examined above lack a complete mathematical statement of the problem and a solving algorithm, which makes it more difficult to clearly understand the issue that is being explored.That is why the promising directions of research are considered to be: -improvement of existing algorithms for solving the problem on plastic behavior of bulk materials; -development of the appropriate software code and its verification.

The aim and objectives of the study
The aim of present work is to improve algorithmic approaches to solving a nonlinear problem on the mechanical behavior of bulk materials by the Drucker-Prager yield criterion.This will make it possible to minimize requirements for computer resources.
To accomplish the set aim, the following tasks have been solved: -to formulate a mathematical model for the elastic-plastic behavior of isotropic bulk material; -to improve a procedure for solving numerically a problem on the mechanical state of bulk material based on the reverse-mapping algorithm; -to compare data obtained by numerical experiments with the data acquired employing commercial software products.

Materials and methods for examining the elastic-plastic state of bulk material
According to the incremental theory of plasticity, a mathematical model of the elastic-plastic behavior of an isotropic bulk material includes the equilibrium equations, a generalized Hooke's law, and geometrical equations [5,6,18 where ij σ are the components of symmetric tensor of stress gain of the 2nd rank, Pa; ρ is the density, kg/m 3 ; i b are the components of gain vector of the mass forces, such as gravity, N/kg; E is the modulus of elasticity during axisymmetric compression, Pa; ν is the Poisson's ratio; δ ij is the Kronecker symbol; 0 ij σ are the components of gain tensor of initial stress, Pa; , ε are the elastic and plastic components of gain tensor of total deformations, respectively; , ij ε i u are the components of displacement gain vector, m.
When employing a criterion of the onset of the Drucker-Prager yield, a condition for the yield of bulk material (plasticity function) is written in the following form [5,6] ( ) where F is a function of the bulk material' surface yield; ( ) ( ) , y c σ φ is the yield limit of bulk material, Pa; c is the adhesion force between granules of bulk material, Pa; ϕ is the angle of internal friction, or the angle of native repose of the bulk material, rad.
If we assume that the surface of fluidity by Drucker-Prager streamlines the surface of fluidity by Mohr-Coulomb [5,6], the expressions for ( ) , y c σ φ and ( ) α φ take the form: Initial conditions for (1), (2): Boundary conditions for (1), (2): -displacement vector gain 0, where S u is the surface (or a point of the surface), on which the displacement is set, m 2 ; -symmetry 0, where n i are the components of external normal vector to the body surface; S su is the surface of body symmetry, m 2 .We shall consider basic theoretical principles of the implicit return-mapping algorithm) [5,6].True elastic stresses in the case of occurrence of elastic-plastic deformations in a bulk material are determined by relation

(
) where σ ij are the components of stress tensor of the 2nd rank, Pa; C ijkl are the components of the fourth rank of material's elastic constants, Pa; tr kl ε are the components of tensor of trial (full) deformations of the 2nd rank, which is determined in the vicinity of elastic medium; ( ) are the components of plastic deformation tensor of the 2nd rank; Δε are the components of gain tensor of plastic deformation at the i-th step of loading; N is the number of loading steps.
For the case of non-associative law of plastic flow at α≠β, we have where ( ) ( ) A scalar associative multiplier Δλ or plasticity coefficient (8) in the absence of strengthening is determined from formula ˆ: : , ˆ: : where 1 ˆ2 1 ˆˆ: 2 are the tensors of the 2nd rank, which are derived from functions ( 2) and ( 9) by a stress tensor, respectively; 4 Ĉ is the tensor of elastic constants of the 4th rank, Pa; Î is the singular tensor of the 2nd rank; ˆtr ε is the tensor of trial elastic deformations at each step of loading.
Considering ( 8) and (10), formula (7) for the loading step k+1 can be rewritten in the form ( ) ˆ: , where 4 ˆ: ε is the tensor of trial stresses that is determined in the vicinity of elastic environment, Pa.
Formula (11) yields a mapping of the trial stress tensor ˆtr σ in the direction of the yield surface.That is why this method of integration, which is built on the inverse Euler method, acquired the name of the return-mapping algorithm [5,6].
The system of equations (11) taking into account the symmetry of stress tensor has 7 unknowns, in particular, 6 independent components ˆtr σ and Δλ.That is why, in order to receive uniqueness, the systems of equations (11) are to be supplemented with the scalar equation ( 2) in the form of a requirement that the condition of yield is satisfied at the end of a loading stage The non-linear system of equations ( 11), ( 12) can be rewritten in the form of discrepancies.In this case, it is required to pass over to a six-dimensional space considering the symmetry of stress and deformation tensors.This makes it possible to replace tensors of the 2nd rank σ n and m with six components.Thus, rather than using a tensor of the 4th rank, one can employ elastic constants tensor of the second rank with dimensionality 6´6: ( ) ( ) , .
In order to solve a system of nonlinear equations ( 13), a Newton's method is typically applied whose iterative procedure is recorded as follows: where ( ) Here, index k refers to the step of loading while index jto the number of iteration by the Newton's method.
At each iteration step by the Newton's method ( 14) it is advisable to not find the inverse matrix but solve a system of linear algebraic equations (SLAE) using the Gauss exclusion method.This makes it possible to significantly reduce the number of arithmetic operations, by about 3n 2 , where n is the dimensionality of SLAE.
For k=1, one solves an ordinary elastic problem relative to total displacements provided the boundary conditions ( 4)-( 6) are assigned, which determine the trial stresses.Next, for the part of the layer of a bulk material that is in the elastically-plastic state, one determines the gains of plastic deformation and the tensor of elastic tension from solution (14) and finds initial stresses from formula The next steps of integration (1), ( 2) for k>1 are performed only with a load by initial stresses (15), (16), that without taking into account external and gravitational loads and at boundary conditions (5).In this case, one solves the elastic problem and determines the gains of displacements Δu k and refines the values of total displacements, which determine the new trial stresses.Then, from solution (14), one determines a new value for the gain of plastic deformation and the tensor of elastic stress for the body part that is in the elastic-plastic state.Next, in order to perform the next step of loading, one finds the gains of initial stresses from formula The criterion for the end of computations can be, for example, fulfillment of condition New trial stresses σ tr (k) in the algorithm of solving the problem can also be determined through the previous values of σ tr(k-1) and the gains of elastic deformations Δε k , which are found via Δu k , from formula In order to determine Δu k , at each step of integration by time, we use a gain of initial stresses in the form Total plastic deformations are determined from formula In the case of the associative law of plastic flow, at α=β (γ=φ), fluidity functions F (2) and G (9) converge.Then m=n, and the direction of gain in plastic deformation at current becomes normal to the yield surface.In this case, it is necessary to replace m with n in formulas ( 6)- (19).In other words, the problem is somewhat simplified.Further algorithm for solving the problem with the associative law of flow is the same as in the case for the non-associative one.
In order to numerically implement the proposed algorithm, we used a finite element method (FEM) and the high-level programming language Fortran employing the integrated development environment Compaq Visual Fortran [19].In this case, a global matrix of SLAE is built in the ribbon form and SLAE is solved by the Gaussian method taking into account the ribbon form.

Results of numerical studies into elastic-plastic state of bulk material
Testing of the developed programming code for solving a problem on the elastic-plasticity of bulk material was performed on the example of a model material, characterized by the associative law of current at different values of the angles of native repose.Construction of the tetrahedron grid was executed in the CAD-system for grid generation Gmsh [20].
Test.A problem on the elastic-plasticity of a bulk material using a classic model of Drucker-Prager.The estimated area is three-dimensional, ¼ of a cone with radius 2 2 0,34 r x y = + = m and height z=0.3 m.Physical properties of the bulk material: apparent density ρ=800 kg/m 3 , modulus of elasticity E=4000 Pa, Poisson's ratio ν=0.45, adhesion force between granules of the bulk material c=400 Pa, angle of native repose ϕ=15, 10, 5°.Load: gravitation g z = =-9.81m/s 2 .The associative law of current α=β.Boundary conditions: fixing on the xOy plane -u z=0 =0, symmetry along the xOz and yOz planes.That is, the terms of the problem, as well as the solution, correspond to a two-dimensional axisymmetric problem.
Results of solving the problem, as well as their comparison with the data obtained using the software ANSYS Mechanical APDL [21] for an axisymmetric geometry, are given in Table 1.
Results of numerical simulation of the problem on plasticity of bulk material with the use of the developed software are shown in Fig. 1.
To visualize the results of calculations, we applied the open-source graphic package for interactive visualization of ParaView [22].An analysis of comparing the results showed that the data from modeling by using the developed software coin-cide with the numerical solutions, which were obtained by employing the commercial software ANSYS Mechanical APDL [21].In this case, the maximal value of error for such magnitudes as u s , σ eqM , el eqM ε and σ m is within 0.25-1.71%, for pl eqM ε is 2.1-5.3 %, which is sufficient enough to perform engineering calculations.
As a result of verification of the modified algorithm for solving a problem (1)- (7), it was established that: -solving a linearized system of equations ( 14) at each step of iterations using the Gauss method instead of determining the inverse matrix makes it possible to reduce the number of arithmetic operations by about 3n 2 for each plastic finite element; -results of solving the problem with initial stresses via absolute values of displacements in line with (15), ( 16) coincide with the results obtained through the gains in displacements according to ( 17)-( 19) under conditions of the test problem.
The shortcomings of the work are, probably, the following: -the study is performed only for a classic model of Drucker-Prager with the associative law of current; -the software is designed only for the application of linear tetrahedron finite elements; -a lack of comparison of the results of numerical simulation with an experiment and the data obtained by using the ABAQUS software package.
The research results are useful when running a numerical analysis of physical fields of electro-thermal equipment, Table 1 Comparison of solutions to the problem on plasticity of bulk material, obtained by employing our own programming code, and using the software ANSYS Mechanical APDL  which uses bulk materials, and are a continuation of previous studies [3,4,18].
Further research may address the ways of solving the problems on plastic behavior of bulk materials in the CAP approximation and extended Drucker-Prager models employing hexagonal finite elements.

Conclusions
1. Based on the classic model of Drucker-Prager, we formulated a mathematical statement of the problem on elastic-plastic behavior of isotropic bulk material.
2. We have improved the procedure for solving numerically a problem on the mechanical state of bulk material using the return-mapping algorithm.The modified technique makes it possible to reduce the number of arithmetic operations by about 3n 2 (n is the dimensionality of SLAE) at each iteration step for each plastic finite element.
3. We have tested the programming code on the example of a model material, characterized by the associative law of current, at different values of the angle of native repose.It was found that the maximal value of error for such magnitudes as u s , σ eqM , el eqM ε and σ m is within 0.25-1.71%, and for pl eqM ε is 2.1-5.3 %.

Introduction
The phenomenon of separation of interlayer bonds is observed in the operation of highways, airfield pavements, foundations of high-rise buildings with supports or isolation joints.The contact between the building sole and the base often leads to emergency situations.Such phenomena are modeled by contact problems with unilateral constraints.
A detailed research of the fundamental laws of contact interaction requires a comprehensive consideration of the geometric features and imperfections of the conjugate surfaces, physical and mechanical phenomena in the areas of their direct attachment (friction, slippage, adhesion, etc.).The numerical solution of contact problems is usually carried out on the basis of the finite element method.For the modeling of unilateral constraints, various physical models are used.The research on the stress-strain state of the model and the state of the contact zone (presence of a friction zone and a separation zone) can be defined as the solution of a direct problem.
Studies of contact problems of mechanics of deformable media are conducted in two directions.Within the framework of the first direction, problems of conjugation of media with sharply differing mechanical properties (layer-inhomogeneous elastic and elastoplastic bodies, conjugation of solids to liquid or gaseous media, etc.) are considered.In such problems, the boundaries of the contact area are, as a rule, specified and not changed in the course of deformation.To solve these problems, the methods of the function theory of a complex variable and the theory of potentials with integral is the equivalent stress by Drucker-Prager, Pa; 1 3 ij ij ij kk s = σ − δ σ are the components of deviatory stress tensor, Pa; m 1 , 3 ij kk σ = δ σ Pa;

γ
is the angle of dilatation, rad.; ˆ.Ĝ ∂ = ∂σ mA potential function of the Drucker-Prager fluidity criterion is expressed by relation[3]

Fig. 1 .
Fig. 1. Results of numerical simulation of the problem on plasticity of bulk material: a -total displacements; b -equivalent plastic deformations by von Mises; c -equivalent stresses by Drucker-Prager; d -hydrostatic pressure u s -resultant displacement, m; σ eqM -equivalent stress by von Mises, Pa; σ eqDP -equivalent stress by Drucker-Prager, Pa; el eqM ε -equivalent elastic deformation by von Mises; pl eqM ε -equivalent plastic deformation by von Mises; σ m -hydrostatic pressure, Pa; ϕ -angle of native repose, degrees; σ y -yield limit of material, Pa