Impact loading of a space nuclear powerplant

Preferred formulation of the problem in two space dimensions are described for solving the three fundamental equations of mechanics (conservation of mass, conservation of momentum, and conservation of energy). Models of the behavior of materials provide the closure to the three fundamentals equations for applications to problems in compressible fluid flow and solid mechanics. Models of fracture and damage are described. A caloric model of the equation of state is proposed to describe thermodynamic properties of solid materials with the phase transitions. Two-dimensional problems of a high-velocity impact of a space nuclear propulsion system reactor are solved. High-velocity impact problems of destruction of reactor are solved for the two cases: 1) at its crash landing on the Earth surface (the impact velocity being up to 400 m/s); 2) at its impact (with velocity up to 16 km/s) with the space debris fragments.


INTRODUCTION
he problem of quantitative description of the impact is a rather complex task, therefore it is connected with a whole scientific and technical research area rapidly developing lately.The bodies collision is accompanied by various processes, which emergence and relative role depend on the shape and physical characteristics of the objects and, more importantly, on the relative velocity of their collision.Moreover, the collision is often accompanied by penetration of one body into the other.In the process of high-speed collision, one or both of the bodies collided can collapse, scatter, or disperse, followed by dissipation of a significant amount of energy.In emergency situations, modern space vehicles with thermionic reactors "shoot off" (jettison) the nuclear powerplants (a simplified sketch of such a powerplant is shown in Fig. 1).There is a certain probability, however, that the reactor fragment containing the nuclear fuel reaches the Earth's surface despite considerable thermal and mechanical loads in dense atmospheric layers.The velocity of the impact of the remaining part of the reactor system can reach 400 m/s.It is next to impossible to solve impact problems of real engineering objects though the computational engineering development level is rather high and fairly realistic mathematical models of material behavior are available.The reasons are the complicated spatial locations of the reactor fragments and the multiscale character of the problem.In such situations, the object considered is simplified, which makes it possible to construct a number of models aimed at studying the influence of the impact parameters on particular basic fragment of the object.The simplification used implies that the materials of small-scale parts were averaged inside the reactor zone in the additive approximation.The mass of the non-principal materials was assumed to be too small (beryllium, uranium dioxide and zirconium hydride compose 95-97 % of the reactor mass) to exert any significant effect on the shock-wave amplitude.As such a medium (mixture) does not have the volume defects, its specific volume on the wave front can be calculated as T mix 1 ( ) ( ), where Vi is the specific volume of the i-th species under shock compression of each species separately, n is the number of species in the mixture, and  i is the mass concentration, m i is the mass of the i-th species.
Thus, our study is based on the assumption that the additivity rule is satisfied rather accurately.In the additive approximation, the volume of the shock-compressed mixture is assumed to be equal to the sum of the volume of the species obtained at the same pressure with their separate shock compression in the form of homogeneous monolithic samples.It was further assumed that the external elements of the structure burn down when the reactor enters the dense atmospheric layers, and the reactor remainder is an object with a complicated internal structure illustrated in Figs. 1.The reactor consists of a beryllium shell, uranium dioxide fuel cells, and zirconium hydride fillers.The end-face (longitudinal) and side impacts were considered.In the first case, we have a problem of an impact of the cylinder side surface onto a deformed target.A specific feature of this formulation of the problem is a multiply connected computational domain with a large number of contact surfaces.In the second case, the reactor model is formed as a ring-shaped structure, while the computational domain is again multiply connected and numerous contact surfaces are formed (Fig. 1).

MATHEMATICAL PROBLEM
he problem was solved in a formulation proposed by Wilkins [1] with the symmetric algorithm of contact boundary calculation [2].Previously, we proposed the process of automatic construction of a triangular grid for an arbitrary multiply connected domain [3][4].According to [1] we used Lagrange coordinates permits the history of mass elements to be followed where the integrated effects of plasticity and external loads change the material physical properties.In order to describe all the dynamic processes, firstly, the main conservation laws must be taken into account (mass, momentum, energy and equation of state), elastic and Prandtl-Reuss stress-strain relations: -equations for tracks of the material particles: ij s : deviator of the stress tensor -P: pressure -: density -u: speed Prandtl-Reuss stress-strain relations 2 , with the plastic strain increments given by Levy-Mises Equation where G -shear modulus, 0 Y -dynamic yield stress, for definition scalar d  known procedure used to bring to the yield circle.The plastic flow will be described by conserving the deviator of stress tensor at the elastic limit.We will assume that after the deformation increase, the stresses change and exceed the yield circle.This contradicts the yield condition and cannot be realized.Instead, we will assume that a plastic flow takes place in the material and the stresses remain at the elastic limit, i.e. on the yield circle.The plastic part of the deformations is perpendicular to the yield surface that leads to limitation of exactly those stresses which are connected with this part of the deformations.Thus, summarizing the assumption about the yield, we can write down the yield condition in the terms of the stress tensor deviator.An algorithm realizing such a condition was proposed by M.L. Wilkins [1].If a redundant change of the stresses in a particular element has led to a law violation, then the main values of the stress tensor deviator have to be corrected in such a way, that the condition is met again.

EQUATION OF STATE
e consider the three-term Mie-Grüneisen equation of state with the solid-phase free energy being determined as x vl ve where V is the specific volume,

 
x E V is the "cold" energy, T is the temperature, , 3 / 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 is the experimental value of the electron heat capacity under standard conditions.The elastic (cold) component of energy Ex(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 parameter where 0 / is the adiabatic modulus of volume compression, v c is the specific heat at constant volume, and ,0 t P is the thermal pressure in the initial state.The general expression for the volume dependence of the Grüneisen parameter has the form In Eq. ( 3), the situation value corresponds to the Landau-Slater theory [5,6] at t=0, to the Dugdale-McDonald theory [7] at t=1, and to the free-volume theory [8] at t=2.
To determine the zero isotherm, we equated the expression for the Grüneisen parameter (2) at the zero temperature ( 0K T  ) to the expression for the generalized Grüneisen parameter (3): Here, x a is the value of the parameter 0 T a  at the zero temperature in Eq. ( 2), which can be taken as the first approximation.The differential Eq. ( 4) has an analytical solution for "cold" pressure and energy: Using the definition of the Grüneisen parameter in the Debye approximation and Eq. ( 2), we obtain the characteristic Debye temperature on the volume: is the Debye temperature at the initial conditions.The constants for Eq. ( 5) were determined and calculated in Fig. 2 and [9].It was also demonstrated there that the set of semi-empirical relations (1)-( 5) 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

MODIFICATION OF THE EQUATION OF STATE.
hermodynamic and kinetic properties of liquids usually cannot be predicted theoretically on the basis of the first principles only.The calculation of phase diagrams is additionally complicated by the fact that the most important characteristics of phase transitions (heat of transition, difference in phase densities, etc.) are small differences in T quantities that have large values and cannot be calculated with sufficient accuracy.Because of thermodynamic nonequilibrium typical for polymorphic transformations in shock waves, the general thermodynamic relations for phase transitions (equality of chemical potential in the region of simultaneous existence of the phases and Clausius-Clapeyron equation) as applied to shock-wave processes can only be used as approximate estimates.Therefore, beginning from the famous van der Waals' paper, numerous attempts have been made to construct the dynamics of a liquid substance by extrapolating the known thermodynamic functions from different areas of the phase diagram.The absence of a commonly accepted thermodynamic model for liquids, which would be equivalent to the Debye approximation for the crystalline state, is a severe obstacle for constructing the equation of state for liquids.In studying particular models of thermodynamic states, it is clear that the usual classification of states in the region of high pressures and temperatures often loses its definiteness and becomes conventional, while the boundaries between the phases either disappear altogether or get fuzzy and actually correspond to continuous mutual transformation of states close to each other.The substance is both compressed and heated in the shock wave.In comparatively weak shock waves propagating over a cold substance, however, the pressure is mainly increased owing to compression.The pressure growth rate in relative units exceeds the temperature growth rate, and the increase in the "cold" compression energy is much greater than the increase in the thermal energy.As the shock wave intensity increases, the relative contribution of the thermal components of pressure and energy increases and becomes prevailing in strong shock waves.To take into account melting, the equation for a solid body (1) is modified.Though the phonon spectrum in the liquid phase is obviously a non-Debye spectrum (as the liquid does not have the far order, its molecules do not have forbidden states; hence, the molecular motion is not discrete), such a modified model can provide positive results at sufficiently high densities of the liquid near the line of phase equilibrium between the liquid and solid phases, because the motion of a liquid molecule within the first coordination sphere at rather high densities can be conventionally considered as vibrational motion.Thus, the higher the pressure, the more precisely the model is satisfied.The expression for the free energy ( , ) L F V T of a monatomic liquid in the so-called rough classical model of the harmonic oscillator has the form where , ( ) x L E V is the "cold" component of the internal energy of the liquid and  is the characteristic "Debye temperature" of the liquid.Formally, the expression for the free energy of the liquid (6) differs from the expression for the free energy of the solid (1) only by the presence of an additional entropy term c F RT A   . Nevertheless, the differences between the solid and liquid states are much more essential.Thus, different cold curves are used to describe the solid and liquid states.First, the major part of the entropy jump is related to the change in the solid structure due to the transition to the liquid state, which is caused by the loss of the far order and leads to the formation of a collective or configuration entropy (equal to zero in an ideal crystal) in thermodynamic models of the liquid.The configuration part of the entropy c c S F T    characterizes the measure of liquid disordering and should remain a finite quantity (similar to the entropy of amorphous solids) as the temperature formally tends to zero 0 T K  .Second, the zero isotherm of the liquid is shifted with respect to the zero isotherm of the solid toward lower densities.The primary reason is that the liquid-phase density extrapolated to the domain of low temperatures is somewhat lower than the density of the solid substance.These principal differences between the thermodynamic descriptions of the liquid and solid states form the basis for the thermodynamic model modification considered in the present work.Thermodynamic functions of both solids and liquids are formally examined in the entire temperature range, including the domain close to 0 T K  .This is convenient for a unified description of both phases and allows functions typical for the solid state (such as the density at the temperature 0K T  , zero isotherm, etc.) to be used to describe the liquid state thermodynamics.We assumed that the physically meaningful branches are those at m T T  for the liquid state and at m T T  for the solid state ( m T is the melting temperature).To "control" the configuration entropy we introduced a certain parameter s a , which has the meaning of the residual entropy at the temperature 0 T K  . Then, the entropy term of the free energy of the liquid in Eq. ( 6) takes the form / , which does not contradict the approximation of the "rough" model of the liquid at 1 s a  .This fitting parameter makes it possible to define correct values of the entropy and specific volume jumps on the melting curve.
Considerations used to derive the expression for the Grüneisen parameter are not confined to the condensed phase.For this reason, all relations in contain only the general thermophysical properties of the material, which are defined and have an identical meaning both for the solid and form the liquid states.Therefore, repeating all transformations applied in solid, we can obtain equations for ( ) x E V , ( ) x P V , etc., whose functional form is similar to those for the solid state.The differences are only in the parameters determined by particular initial conditions.Therefore, the use of the modified equation of state requires only six constants, as the equation for the solid phase.The constants for several liquid metals can be found in.Moreover, the liquid model requires the jumps of the volume V  and entropy S  due to melting to be known.

MELTING AT HIGH PRESSURES BEHIND THE FRONT OF A STRONG SHOCKWAVE
he phase transition of the first kind (melting) is understood as an equilibrium transition of the substance from one phase to another with jump-like changes in the first derivatives of the Gibbs energy G with respect to temperature and pressure, i.e., the entropy S and specific volume V experience jump-like changes during melting: The entropy of the liquid phase is always greater than the entropy of the solid; therefore, the entropy change during melting is always positive.The change in the volume during melting, however, can be either positive or negative.The inequalities 0 S   and 0 V   , which mean greater orderliness and density of the crystalline phase as compared with the melt, seem to be natural.They are valid for most substances and ensure a positive slope of the melting curve 0 dP dT  .At the same time, there are some substances (e.g., gallium, bismuth, and water) with negative values of this derivative 0 dP dT  .As the shock-wave pressure increases, the thermal energy imparted to the substance continuously increases and the transition of the initially solid substance to the liquid state is expected to start at a certain pressure level.The further behavior of the dependence ( ) T P along the dynamic adiabat can be understood by analogy with melting at atmospheric pressure, where an increase in the energy imparted to the substance starting to melt does not lead to an increase in temperature until the substance becomes completely melted.Further heating of the liquid is accompanied by an increase in temperature.A similar pattern should also be observed under shock compression, with the only difference that a certain increase in temperature can be expected in the domain of simultaneous existence of two phases (segment of the melting curve between the shock adiabats for the solid and liquid substances, in Fig. 3), because / 0 dT dP  for the melting curve of the majority of substances (so-called "standard" substances).
The last equation in (7) is the equality of the chemical potentials of both phases per 1 mole of the substance.The subscripts refer to the solid state (S), liquid state (L), and values on the melting curve (m).

CALCULATION RESULTS OF EOS
sing the system of equations of state of the solid (1) and liquid (6) phases and taking into account the phase equilibrium condition (7), we calculated the dynamic adiabats and the melting curves for several metals.The numerical calculations were performed by an iterative method.The shock adiabat of the solid body was used as the reference curve.In this case, we first assumed that 1 s a  in Eq. ( 6) and chose a value of x a in Eq. ( 4) that ensured the minimum difference between the calculated and the reference shock adiabats of the liquid.After that, for this value of x a , we found the value of as from the condition of identical chemical potentials of the solid and liquid phases for specified entropy and volume jumps on the melting curve.After that, the procedure was repeated with the new value of s a .The calculated adiabats with the phase transition in the T P  coordinates are shown in Fig. 3.The correctness of taking into account melting behind the shock-wave front in the equation of state derived is indirectly confirmed in the paper of Sakharov [11] who measured the viscosity behind the shock-wave front in aluminum and concluded that it remains in the solid phase up to pressures of 1 Mbar.The calculations show (see Fig. 4) that aluminum melting begins at the pressure P=1.07 Mbar.The transition of Al to the liquid state is finalized at the pressure P=1.2 Mbar.The calculation for lead shows (see Fig. 4) that Pb melting begins at the pressure P=0.36 Mbar.The transition of Pb to the liquid state ends at the pressure P=0.41 Mbar.The result obtained is in reasonable agreement with the data [12], where it was noted that melting in the shock wave begins when the mass velocity reaches ~ 650-700 m/s, which corresponds to pressures of 0.23-0.25 Mbar.In [13] experiment, lead behind the shock-wave front is already in the melted state at pressures of 0.4 Mbar. Figure 5: Shock adiabats of metals with allowance for melting.
In accordance with Urlin theoretical predictions [14], we performed calculations (Fig. 5), which show that melting exerts a minor effect on the adiabat shape in the P V  coordinates.The inflections on the diagram in Fig. 5 are so small (the melting zone is indicated by the asterisks) that they can be hardly distinguished visually.The values of these inflections are comparable with the error of data obtained in shock experiments.Now let us consider the melting curve found as the boundary between the phases with the corresponding equations of state.Fig. 6 shows the melting curve calculated for aluminum, the experimental data and Simon's melting curve.The calculated melting curve lies somewhat lower than the optimal curve, and the error of calculating the temperature on the melting curve is smaller than 10% in the entire range of pressures considered.U Fig. 7 shows the lead melting curve and the experimental data; the results are seen to be in good agreement.As an estimate, the figure also contains Simon's melting curve for lead extrapolated to the range of high pressures.As the parameters were chosen on the basis of low-temperature experiments, Simon's melting curve underestimates the melting temperatures in the high-pressure domain.
Comparing the calculated results with available experimental data, we can conclude that the calculation accuracy is fairly high and the phase equilibrium condition (7) can be used for calculating the melting curves for the metals considered.

FRACTURE
inematic strength characteristics include the limiting values of elongation (usually under uniaxial tension) and shear.Brittle materials are also destroyed by compressive strains.Kinematic characteristics are accumulated quantities incorporating the entire history of the process.The most frequently used quantities are the limiting elongation and shear strain.Finding these quantities involves calculation of the primary tensile and compressive strains and also the shear strain If the tensile strains in the course of deformation exceed the limiting elongation * 1

  
or the limiting shear strain * ,   then the element material is assumed to be fractured, i.e., it has no longer resistance to tension and shear, but still has resistance to compression.Force strength parameters include the limiting values of tensile [20], compressive, and shear stresses.If the stresses proper are used, then they are instantaneous criteria, i.e., as soon as the principal stresses exceed the limiting values, the element material is assumed to be fractured: Most material, however, possess properties of plasticity and viscosity; therefore, their fracture requires a time interval during which the material is under overstrain.The code involves one of such criteria (it is demonstrated by an example of the principal tensile stress) [21]:

CONVERSION OF FRACTURED ELEMENTS TO PARTICLES
f the element is at the boundary of the computational domain and the parameters of the material reaches a limited value, the material of the element is replaced by discrete particles.Radius of the particle is calculated from the condition of incorporating one or more of the particles in the element.Mass of the element is allocated between the discrete particles.For one time step is only one layer of boundary elements can be converted into discrete particles, as is that the velocity of the wave front of destruction does not exceed the speed of disturbances in the medium [2].Fig. 8 shows conversion of triangle elements (A, B and others) to particles [22]:  Element A is removed from the element grid;  Particle A is added as a particle node;  All of the element variables (stress, strain, damage, etc.) are transferred to the particle;  The mass, velocity and center of gravity of the particle node are set to those of the replaced element.The nodal velocity is obtained from the momentum of the element (three nodal masses and velocities);  The masses of nodes b, c and k are reduced by the removal of element A;  For the conversion of element B (which has two sides on the surface) to node B; most of the steps are similar to those used for element A.

TWO-DIMENSIONAL PROBLEM OF THE IMPACT OF A MODEL NUCLEAR POWERPLANT "TOPAZ" ONTO THE EARTH'S SURFACE
s the first example, we consider an impact of the reactor on sandstone with the initial velocity of 400 m/s in the axial formulation.Realistic photo of the nuclear powerplant with thermionic reactors for space applications "Topaz" is illustrated in Fig. 9. Geometry of the reactors is illustrated in Fig. 1.In the case of the reactor impact I A onto sandstone, the wave pattern is complicated by the fact that sandstone is destroyed by compressive stresses.This specific feature is induced by the inner structure of sandstone where strong crystals of sand are bound with brittle cement mass.As sand and cement have different compressibility, the shock wave forms shear stresses on interfaces between the media, which destroy the connections on the boundary, i.e., the resultant product is sand with a fine fraction of cement.
Free sand exerts practically no resistance to shear strains.Thus, the impact compression forms a domain of a fractured material near the reactor-sandstone contact surface and, as a consequence, an unloading wave.The interaction of the side unloading waves and the unloading wave from the fracture zone leads to formation of a zone of tensile stresses with a higher amplitude, leading to fracture of the reactor materials (zirconium hydride filler and fuel cells).Beryllium shell fracture follows the mechanism of shear-induced quasi-static fracture.Fig. 10 shows the frames of the impact of the reactor model onto sandstone.Let us further consider the computation of the side impact of the nuclear powerplant.Fig. 11 shows the frames of the computed side impact of the reactor on a granite plate with an impact velocity of 400 m/s.This computation shows that the side impact of the reactor is more critical in terms of the impact velocity than the end-face impact because of higher strains in the contact area.While moving, the reactor filler (which is heavier) weighs on the beryllium shell, thus, inducing shear fracture in the contact area with the granite plate and inner fracture of the casing material, which violates the reactor proofness.Numerous cracks appear in the ZrH filler, which results in disintegration of fuel cells and, as a consequence, possible radioactive contamination of the place of incidence of the reactor.
It should be noted that the results of calculations in the 2D planar case show a loss of symmetry in space.This is because the build perfectly symmetrical grid for complex bodies is not possible.Enough small deviations from symmetry lead to small differences in the strains, but near the critical values of this causes no symmetry of fracture in the space.To impact at the normal calculation can be performed in one half of that mathematically ensures symmetry of fracture.However, in practice we generally are not spatially symmetric picture of damage.

CALCULATIONS OF REACTOR DESTRUCTION DUE TO ITS COLLISION WITH SPACE DEBRIS ON THE ORBIT ("TOPAZ")
n the first calculation, a two-dimensional reactor model interacts with a steel object 3 cm in diameter with a velocity of 11.7 km/s.This velocity is the most probable one for collisions with space debris fragments.The calculation is performed for a reactor model protected by the Whipple shield (aluminum sheet 2mm thick surrounding the reactor).The experiment shows that the kinetic energy of the object with such a velocity is sufficient for penetration through the target, acceleration of the expelled mass, and fragmentation of the impinging object and expelled part of the target.Shock waves arising during target penetration are reflected from free surfaces as unloading waves whose interference generates strong tensile waves in the materials of the Whipple shield and the impinging object.The parameters of the resultant waves are substantially higher than the strength characteristics of the materials, which leads to fragmentation of the structure and to formation of a cloud of fragments behind the target.The further interaction between the cloud and the reactor model proceeds in accordance with the following scenario: the cloud of fragments of the destroyed projectile almost with the initial impact velocity reaches the beryllium shell of the reactor.The remaining part of the cloud also contacts the reactor shell with a certain delay.Despite the presence of the Whipple shield, the kinetic energy of the cloud I of particles is sufficient for destroying the reactor shell made of beryllium.Multiple interactions of particles form a stable spherical shock wave, which destroys not only the shell, but also the zirconium hydride interior and uranium fuel cells.Fig. 12 shows the frames of the interaction of the cloud of fragments with the reactor.The presence of beryllium cylinders in the shell structure is manifested only as an additional stress concentrator and results in more expressed destruction on the shell edges, but exerts practically no effect on the overall pattern of shell failure.After that, the interaction of a steel object 3 cm in diameter impinging along the reactor axis with a velocity of 11.7 km/s is calculated.Fig. 13 shows the photographic records of the impact process.A powerful spherical shock wave propagates from the contact surface.Free surfaces at the periphery of interacting bodies allow the materials of these bodies to unload by acquiring velocities directed away from the "center of pressure", thus forming zones of tensile strains and stresses.The destruction process begins.As the size of the impinging particle are small, as compared with the reactor size, an analogy can be found with a point source of pressure applied on the boundary of interaction of bodies, moving along the reactor axis.As the central part of the reactor has a complicated structure, multiple interactions of compression-unloading waves lead to its failure even at the stage of overall compression behind the shock wave front (mainly due to shear strains).The acquired velocity induces systematic tensile strains, which rapidly reach critical values, and fragmentation of the materials of the central part of the reactor.Massive deformation of the central part appreciably alleviates peripheral loads, which allows the beryllium shell to remain in the non-destroyed state.

CONCLUSIONS
1. Calculations show that the high-speed impact of space debris leads to full failure of the reactor.2. Impact of the reactor onto Earth's surface at a speed of 400 m/s leads to decompression of the reactor and radioactive contamination of the fall.
values of these quantities under standard conditions, which can be found in reference books on physical and mechanical properties of substances.

Figure 3 :
Figure 3: Melting curve and shock adiabats in the domains of the solid and liquid phases.

Figure 4 :
Figure 4: Shock adiabats of metals with allowance for melting.Figure5: Shock adiabats of metals with allowance for melting.

Figure 6 :
Figure 6: Melting curve for Al.The solid curve is the present calculation; the experimental data  and  [15, 16].The dashed curve is Simon's melting curve.

Figure 7 :
Figure 7: Melting curve for Pb.The solid curve is the present calculation; the experimental data are indicated by  ,  , and  [17-19].The dashed curve is Simon's melting curve.

Figure 8 :
Figure 8: Conversion of triangle elements to particles: (a) interface nodes and elements before conversion; (b) interface after conversion of elements.

Figure 10 :
Figure 10: Impact of a reactor model with a sandstone plate surface.

Figure 11 :
Figure 11: Frames of the calculated side impact of the reactor model onto the granite plate surface

Figure 12 :
Figure 12: Photographic records of the impact of space debris with a reactor equipped with the Whipple shield.

Figure 13 :
Figure 13: Photographic records of the impact of space debris with a reactor (axial formulation).