Numerical simulation of deformation and failure processes of a complex technical object under impact loading

The main points of development of numerical tools for simulation of deformation and failure of complex technical objects under nonstationary conditions of extreme loading are presented. The possibility of extending the dynamic method for construction of difference grids to the 3D case is shown. A 3D realization of discrete-continuum approach to the deformation and failure of complex technical objects is carried out. The efficiency of the existing software package for 3D modelling is shown.


Introduction
Numerical modeling of the matter behavior under the influence of shock load requires taking into account the laws of elastoplastic flow and, therefore, of stress and strain. Such problems are solved using the equations of mass, momentum, and energy conservation [1][2][3][4]. The physics of the processes at high-speed collision of solid bodies can be studied in one-dimensional case, for example, construction of equations describing the state of materials. However, there is a need for both two-dimensional and spatial calculations when modeling the processes of deformation and destruction of complex technical objects that do not fit into one-dimensional schemes. Thus, the two-dimensional numerical tools cannot be used to describe the phenomena of a nutating cylinder rebound off an obstacle in collision at an angle or the impact of a complex technical object, for example, a model reactor at its center line due to its multiply connected structure. Therefore, there is an urgent need for further development of the existing two-dimensional numerical tools [5] to be used in the spatial case.
The software package is designed to solve the time-dependent problems of deformation and fracture of materials at loading of complex technical objects. The package implements a continuous-discrete approach, combining Lagrangian coordinates for the solid material and discrete particles of finite size simulating fragments of the fractured material. The behavior of the materials is described by an elastoplastic model with a few-parameter equation of state. The difference ratios are put on a tetrahedral grid, which is dynamically constructed in complex multiply connected computational domains. All this allows one to simulate the processes of impact loading of complex technical objects to "the end," i.e., till the complete absorption of the kinetic energy of the projectile body and/or splitting of the object into separate fragments.
The software package developed for the spatial modeling includes: • equations of mass, momentum, and energy balance [6,7]; • the initial data and boundary conditions for all elements of the object [3,6,7]; • equations of state and equations of the process for the reactor materials [8][9][10][11][12]; • necessary criteria for the destruction of materials covering a variety of mechanisms of their destruction [2,3,5,13]; • algorithm for replacing cells containing the fractured material by discrete particles of finite size [3,5,7]; • extension methods for dynamic construction of the difference grid [5,14] to the spatial case; • finite-difference equations for the balance equations on a tetrahedral grid [7]; • symmetric algorithm for calculating the contact boundaries and condition for obtaining stable solutions of the system of difference equations [5,7].
Let us consider the above steps in more detail.

Formulation of the problem
We start from the Lagrangian method for describing the environment as the most appropriate way to describe the interaction of deformable solids [1]. The model consists of equations of motion, mass, momentum, and energy balance, as well as the equations of state and elastoplastic Prandtl-Reuss flow: equations of the trajectory of motion of material particleṡ equation of the medium continuity law of change in the material particle momentum change of the internal energy of a particle the strain rate tensor has the formε the stress tensor can be presented in standard form where s ij is deviator of the stress tensor responsible for the reaction to the shear deformation of the material particles; δ ij is the Kronecker delta; P is a function of pressure in Mie-Grüneisen form.
The process equations are taken in the Prandtl-Reuss form where Y 2 0 is the dynamic yield strength, and the scalar factor dλ ′ is determined by the known procedure of reducing to the yield circle. The above equations use the standard notation: each of the indices i, j takes values of 1 and 2; the repeated indices assume summation; the dot over a symbol is the time-derivative; the index after comma is the derivative with respect to the corresponding coordinate; x i , u i are components of the vectors of position and velocity of a particle, respectively; ρ is the current density; G is the shear modulus.
The initial conditions for the ith body at t = 0 in the domain D(x, t) often have the form where ρ i0 (x), u i0 j (x) are specified initial distributions of the material density and the velocity vector throughout the region D(x, t).

Equation of state
We consider a three-term equation of state with the solid-phase free energy determined as where V is the specific volume, E x (V ) is the "cold" energy, T is the temperature, c v,l = 3R/A is the specific heat of the lattice at constant volume, A is the mean atomic weight, R is the gas constant, θ(V ) is the Debye temperature, and c v,e0 is the experimental value of the electron heat capacity under standard conditions. The elastic (cold) component of energy E x (V ) is related exclusively to interaction forces between the body atoms and is equal (including the energy of zero vibrations) to the specific internal energy at the absolute zero temperature. The thermodynamic model of a few-parameter equation of state is based on the dependence of the Grüneisen coefficient on the volume where γ s = βK s V 0 /c v , K s is the adiabatic modulus of volume compression, c v is the specific heat at constant volume, and P t,0 is the thermal pressure in the initial state. The general expression for the volume dependence of the Grüneisen coefficient has the form In equation (12), the situation value corresponds to the Landau-Slater theory [15,16] at t = 0, to the Dugdale-McDonald theory [17] at t = 1, and to the free-volume theory [18] at t = 2.
To determine the zero isotherm, we equate the expression for the Grüneisen coefficient (11) at zero temperature (T = 0K) to the expression for the generalized Grüneisen coefficient (12) Here, a x is the value of the parameter a T =0 at the zero temperature in equation (11), which can be a x = a(0) = 1 + 2/(γ s − 2/3) as the first approximation. Differential equation (13) has an analytical solution for "cold" pressure and energy: where C 1 , C 2 , and C 3 are constants of integration and H 1 (V ), H 2 (V ) are polynomials computed by the formulas [10] The following replacements are made here: We use the definition of the Gruneisen coefficient in the Debye approximation γ = −(d ln θ/d ln V ) T and equation (11) to obtain the characteristic Debye temperature dependence on the volume: where L 0 = L(V 0 ) is the Debye temperature under the initial conditions. We also demonstrated there that the set of semi-empirical relations (10)-(13) describes the behavior of thermodynamic properties of solids within 5-10% in a wide range of pressures and temperatures. For the equation of state to be applied, it is sufficient to know only six constants V 0 , β, K t , c v , L 0 , and c v,e0 corresponding to the values of these quantities under standard conditions, which can be found in reference books on physical and mechanical properties of substances.
The above-constructed equation of state was used to construct the principal Hugoniot for various materials. Figure 1 shows the principal Hugoniot and zero isotherms of aluminum, copper, and lead calculated by the method proposed by the author. For comparison, the results of experiments of various research groups are shown, whose data are put together in [19], as well as the results of calculation of the "cold" pressure [20] (these calculations are classical in the field of creating the equation of state). As one can see, the differences between the "cold" curves calculated in [20] and the authors' dependence are insignificant up to the compression ratios V 0 /V = 1.5. At the compression ratios V 0 /V > 1.5, the differences are more significant (but they do not exceed 7%) due to the method used to construct the model, where the reference state for the construction of equations is the normal state of the condensed matter.
The calculations were carried out for the parameter value t = 0 in equation (16), i.e., according to the theory of Landau-Slater. The calculations for the uranium dioxide UO 2 showed that there are no values of the fitting parameter and/or the initial data, which determine the initial course of the principal Hugoniot, that would allow approximating the experimental data from [20] at t = 0. For t > 1, the approximation is possible, but the derivative ∂P/∂ρ at the starting point for these experiments is quite small. Figure 2 shows the calculation of the shock Hugoniots at t = 0 by the Landau-Slater method and at t = 0.6 in equation (16). Figures 3 and 4 show how the few-parameter equation of state is successfully used to calculate the behavior of ceramic materials (B 4 C and Al 2 O 3 ) at high pressures and temperatures.

Verification of the tool for high-velocity interaction and damage of solids
A test simulation was made to verify the program. The initial conditions and experiment results [34] are shown in figure 5. The experiments were carried out with tungsten alloy rods of equal length-to-diameter ratio L/D = 20 and equal mass m = 74.2 g at a nominal impact velocity of 1.615 km/s for the cylindrical rod penetrator with angles α 1 /α 2 = +2.4 • / − 0.7 • .
An example of conversion computations is shown in figure 6, where a tungsten rod impacts and perforates a steel plate. There is a good general agreement between the computed and 6 1234567890 ''""   experimental results. The nose of the rod is eroded, the fragment distribution is similar (on both sides of the plate), and the hole in the plate is similar. The failed elements are replaced with discrete particles during the computation. If an element is at the boundary of the computational domain and the parameters of the deformation reach a certain limit value, then the element is replaced by a discrete spherical particle. The radius of the particle is calculated from the condition of lying inside the element. The mass of the element is transferred to the particle. Only one layer of boundary elements can be converted into discrete particles in one time step, so that the velocity of the destruction wave front does  not exceed the speed of disturbances in the medium. The conversion of triangle elements to particles is presented [3,7]. For visualization purposes, the particles from the tungsten alloy rod are shown in red in figure 6, and the particles from the steel plate are shown in blue. When the rod hits the plate at an angle, the deformation process starts with formation of a cylindrical track ahead of the front edge of the rod. A "wave" of displaced plate material is formed ahead of the blunt tip of the rod, which then splashes out when the tip plunges into the target (figure 6, 40 µs time).
At the same time, the stresses increase and spread as a wave from the contact surface and transfer the momentum to the target material. When the compression wave hits the backside of the plate, the unloading wave is formed and the speed of the plate material increases in the normal direction. This normalizes the motion of the impacting rod and bends it. This effect is weak when the plate is thin and can only be seen near the tip of the rod in figure 7. The rod makes a through hole in the plate and forms a cloud of fragments. The simulation results are compared with the experimental data in figure 7 and a good qualitative and quantitative agreement is shown.

Examples of application of the tool for high-velocity interaction and damage of solids
To illustrate the efficiency of the numerical 3D tools, we consider the interaction of a model reactor of thermionic power plant with 12 regulating cylinders in the side reflector of neutrons and concentric rows of power-generating channels in the moderator of zirconium hydride. The cylindrical system is formed by fuel elements with integrated thermionic converter of thermal energy of nuclear reaction to electric energy. A great many contact surfaces and various materials show that the software package allows one to solve the problem "completely," i.e., till the complete failure and shutdown of the model reactor fragments. Figure 8 is an illustration of the geometric model of the reactor.
For a heavy tungsten rod, it is easy to penetrate a light and brittle beryllium cap (figure 9, 25 µs). The rod has gone through the cap and its wearing out has barely started. However, when it contacts with heavier cylindrical systems, it quickly loses mass and its length is reduced to a third of its initial value (figure 9, 65 µs time). There are many contacts between different parts of the reactor, which give rise to a complicated pattern of wave interaction consisting of loading and unloading zones. Both the materials of zirconium hydride and homogenized cylindrical systems have low strength which makes the failed material zone spread in the direction of the rod motion, which can be seen in figure 9 (from 25 µs to 100 µs). It should be noted that the impact of the rod on the reactor leads to the catastrophic outcome for the latter.

Conclusions
It is shown that the existing software package REACTOR3D has evolved and is suitable to solve the impact problems in the 3D case. A method for dynamic construction of a tetrahedral mesh for complex, multiply connected objects is presented. An algorithm for modeling a material fractured by discrete particles of finite size is exposed. The key moments of interaction of the fractured material fragments with the boundaries of a solid material and with each other