A computational model for assessing high-velocity debris impact in space applications

Man-made space debris is dominating the background meteorite environment with a growing debris population leading to increased collision risks for satellites, especially in the low Earth orbit and geostationary orbit protected environments. Here we present a computational model for estimating the effect of hypervelocity impact from debris particles on non-shielded propellant and pressurant tanks. Eulerian hydrocode simulation is utilised to model firstly penetration and shock wave formation in the propellant and secondly subsequent detonation wave propagation and interaction with the tank wall. Furthermore, reactive molecular dynamics is used to estimate the risk of detonation in a liquid hydrazine layer. We present simulations of a 3.5 mm aluminium spherical debris particle at a velocity of 14 km/s relative to a hydrazine tank. We find that the degree of damage is strongly dependent on tank temperature and hence on the satellite thermal configuration at its end of life.


Introduction
To mitigate the generation of debris in space, the European Space Agency (ESA) adopted in 2012 a set of requirements applicable to spacecraft operated in low Earth orbit (LEO) Communicated  and geostationary Earth orbit (GEO). As a consequence, at their end of life, spacecraft must be passivated from all stored energy (batteries, propulsion system, reaction wheels) and placed in a safe state. In particular, propulsion systems shall be depleted or made safe, which implies that specific operations need to be planned to deplete the pressurised tanks as far as technically feasible. A combination of thermal drift and collision with man-made debris is the main risk of debris generation identified for a non-passivated propulsion system after it has been decommissioned.
Here we present a physics-and chemistry-based computational model for assessing the risks and potential consequences of such an event. A worst-case scenario has been identified for a spacecraft in the most debris-populated LEO orbit at 800 km altitude and 98.1 • inclination.
The paper is organised as follows. In Sect. 2, the computational methods are briefly described. Section 3 reports the simulation set-up, followed by the main section where the results are presented and discussed. In Sect. 5, we summarise our main findings with respect to tank damage. Finally, in Sect. 6, conclusions are presented.

Finite element models
In this work, LS-DYNA and the FOI-developed code GRALE2D are used for the explicit finite element method (FEM) simulations. The latter includes the SESAME equation of state (EOS) package as an option. The term hydrocode is used interchangeably for these programs. For the debris impact and subsequent plume formation, the forces induced in the material are much larger than the material strength (yield strength) and the deformation is thus dominated by the  [1] material's equation of state. Figure 1 shows a comparison of the two types of equations of state used in this study with data from shock wave experiments [1]. The Mie-Gruneisen equation of state for aluminium, which is used in most of the simulations, shows excellent agreement with experimental shock wave data up to an increase in density of about a factor of two (corresponding to a pressure of a few 100 GPa). A close-up of the region in Fig. 1 where the two models diverge is shown on a linear-linear scale. The SESAME equation of state, which consists of a tabulated data library based on quantum mechanical calculations, agrees well with the experimental data for even higher compression. Furthermore, the SESAME equation of state includes phase transformations and can be expected to provide a more accurate estimate of the material temperature. Figure 2 shows the geometry of a typical impact simulation at different times. Cylindrical symmetry is enforced with the axis of symmetry along the path of the particle. We use an Euler formulation where the material is advected through a stationary or uniformly moving mesh using a second-order advection scheme. The size of the computational elements is 0.2 mm in the vicinity of the impact region. The top left subfigure shows the spherical debris particle (light grey) just before impact against the tank wall (dark grey). Behind the wall is a layer of liquid hydrazine (magenta) possibly followed by hydrazine vapour. This liquid layer can be expected to be thin at end of life. The top right subfigure shows a closeup of the set-up with the computational mesh visible in the background. The bottom subfigures show snapshots during penetration (left) and plume expansion (right) at times 1 and 5 µs after impact, respectively.
For slower processes such as tank wall deformation after penetration and impact at lower velocities, the material model will be important. In this study, we use the Johnson-Cook which couples plastic strain hardening, strain rate sensitivity and temperature softening. A, B, n, C, and M are constants obtained through experiment. The tank material in this study is alpha-annealed titanium alloy Ti-3Al-2.5V. The Johnson-Cook parameters are taken from Zhang [2] where fitting is performed against data from published material characterisation [3]. Equation of state data for liquid hydrazine was taken from Garcia and Persson [4].

Molecular Dynamics
Chemistry at extreme conditions such as high density, temperature, and pressure in the GPa range is very challenging to describe at the atomic level. A relatively new feature in the computational method molecular dynamics (MD) is the possibility to include the breaking and formation of covalent bonds. In this study, we use the ReaxFF [5] force field and the MD-solver LAMMPS [6] to follow shock compression and decomposition of liquid hydrazine. The force field was initially developed for nitramine explosives [7], but has been successfully applied to hydrazine by Zhang et al. [8]. A detailed validation of this type of force field has not been experimentally feasible yet. Nonetheless, predictions of chemical kinetics by ReaxFF have shown good agreement with DFT-based (density functional theory) calculations [5,7]. In this study, we utilise the multi-scale shock technique (MSST) [9] to simulate shock compression by updating the cell volume and temperature in such a way as to restrain the system to the shock Hugoniot and the Rayleigh line. Prior to the MSST simulation, the system was equilibrated using NPT simulation (enforcing constant pressure and temperature) with the Nose-Hoover thermostat/barostat [10]. During the equilibrium phase, the ReaxFF-lg correction [11] was used to calibrate the long-range component of the force field (London dispersion). In this formulation, the correction term has low gradients (lg) at valence distances leaving the already optimised valence interactions intact, but behaves as 1/R 6 for large distances. After calibration, the simulated liquid density agreed with the experimental value of 1.021 g/cm 3 to within 1%. Only two parameters were altered: the dispersion energy correction parameters (see [11]) for the N-N and H-H interactions (set to 400 and 40 kcal/mol A 6 , respectively).

Thermochemical calculations
The thermochemical program Cheetah 2.0 [12] is used to calculate detonation properties of hydrazine vapour. The code solves thermodynamic equations between product species to find chemical equilibrium. Hydrazine data can be found in the standard reactant library.

Set-up
Using ESA MASTER 2009 tool [13], the probability of impact with a standard tank over 25 years (maximum re-entry time) is calculated to be: -10 −2 for a 3.5-mm-diameter aluminium debris at 14 km/s relative velocity. -10 −3 for a 13-mm-diameter aluminium debris at 14 km/s relative velocity.
We consider a typical hydrazine (N 2 H 4 ) monopropellant propulsion system with a 50-cm-diameter titanium diaphragm tank and typical end of life residuals stored at 5 bar and 20 • C. A drawing of a typical LEO tank at beginning of life (left) and end of life (right) is shown in Fig. 3. It is assumed that the tank is externally mounted and therefore does not benefit from spacecraft shielding. A thermal worst-case analysis shows that a maximum temperature of 170 • C can be expected at the tank when the radiator side of the spacecraft is oriented towards the sun and with strongly degraded thermo-optical properties of the materials. In an attempt to justify that spacecraft as passivated today is in a safe configuration, the effect of 3.5 mm debris penetrating the tank and potentially triggering detonation in the residual hydrazine layers is modelled and analysed. In addition to the 3.5 mm debris particle, two other cases are included: -13 mm debris particle at 14 km/s to investigate the effect of increased/prolonged pressure in the liquid hydrazine layer. -0.5 mm debris particle at 5 km/s to establish an approximate threshold for penetration.

Impact and early deformation dynamics
The first step in the analysis deals with tank wall penetration and deformation of the debris particle as it impacts and passes through the wall. These simulations are performed using the two hydrocodes LS-DYNA and GRALE2D. Figure 4 shows a snapshot of the 13 mm debris particle 1.4 µs after impact. The particle travels upward in the image, and the tank is visible in light orange. It has expanded laterally, and a tail of titanium alloy (tank) and aluminium (debris) has formed. Upon impact, the pressure in debris particle and tank wall reaches a peak value of around 3 Mbar corresponding to a density of 6.6 g/cm 3 . During penetration, the pressure wave traversing the particle has a peak pressure of about 1.5 Mbar at a density of 4.8 g/cm 3 . These values are consistent with the comparison presented in Fig. 1 and in line with simulations performed earlier by Kamoulakos et al. [14]. From the SESAME equation of state, we can follow the temperature field in the debris particle during penetration. From the simulation visualised in Fig. 4, we see that the temperature varies between 1000 and 5000 K, implying that the material consists of a combination of molten and evaporated aluminium.  The simulations are performed with the Mie-Gruneisen and SESAME equations of state Such a mixed liquid and vapour state will affect the ejecta dynamics [15]. In order to estimate this effect at the front of the debris plume, simulations with both equations of state were performed. Peak pressure in the vapour and velocity at the plume front were compared 4 µs after impact and are presented in Table 1. Both pressure and velocity are about 20% lower in the SEAME EOS simulation. Even though material melting and evaporation can have a significant effect on ejecta dynamics in general, these simulations suggest that it has a limited effect on the pressure wave that is generated in the hydrazine vapour during the first few microseconds after impact at 14 km/s debris velocity. We conclude that after penetration of the tank wall, the particles will consist of material in the fluid state which can be expected to expand into a plume of material that will interact with the contents of the tank and possibly with the rear tank wall. A study by Kamoulakos et al. [14] indicates that effects such as spalling are small under these conditions and that crack growth after penetration may be suppressed.
Impact simulation of the 0.5 mm debris particle at 5 km/s indicates no penetration of the tank wall, but results in a deep (about half thickness) crater.

Liquid hydrazine initiation
The violent impact and subsequent expansion of material behind the front tank wall will lead to the formation of strong shock waves reflected and transmitted at the material boundaries. For initial tank temperatures below the propellant boiling point, a layer of liquid hydrazine can be expected behind the tank wall. In this section, we investigate the possibility of initiation in such a liquid layer.
There are experiments and calculations suggesting that liquid hydrazine can detonate at a shock wave pressure higher than what can be generated with conventional explosives [16]. Hypervelocity debris impacts can result in such extreme pressures. Figure 5 shows the simulated pressure in a liquid hydrazine layer during penetration of the tank wall. Highpressure regions (in red) can be identified besides the central part of the deforming debris particle. In this region, the pressure reaches values of 10-30 GPa. The region is limited in extension (fraction of mm) and lifetime (few 100 ns). It can be noted that this pressure corresponds to the detonation pressure measured in the reaction zone of conventional explosives which typically can be initiated at pressures in the range of 1-10 GPa.
Reactive molecular dynamics (RMD) was used in order to investigate the dynamic and kinetic response of liquid hydrazine in a microscopic element subject to shock compression in the high-pressure region in Fig. 5. Figure 6 shows visualisations of the hydrazine liquid structure before (top) and after (bottom) shock compression that corresponds to a shock wave velocity of 8 km/s propagating in the direction indicated by the blue arrow. The duration of the transient shock process is about 1 picosecond. Figure 7 shows the pressure (top) and temperature (bottom) over 40 ps after compression. The computed pressure indicates that a shock wave speed of about 8 km/s results in a pressure that corresponds to that seen in the hydrocode simulation. The bottom plot in Fig. 7 shows that the ini- White and green colours correspond to nitrogen and hydrogen atoms, respectively tial compression results in an increase in temperature from 300 K to about 1400 K with no evidence of strong exothermic reactions. This can be confirmed by Fig. 8 (top, blue) which shows that only about 25% of the hydrazine molecules have decomposed 40 ps after shock compression. A shock wave speed of 6 km/s results in a pressure of 15 GPa and a temperature of 800 K. Figure 8 (top, green) shows that this state yields insignificant decomposition. A shock wave speed of 10 km/s corresponds to a pressure of about 47 GPa and a temperature of 2300 K after compression. In contrast to the simulations with lower shock wave speed, the temperature starts to increase after about 5-10 ps. Figure 8 confirms that considerable amounts of NH 3 (about 500 molecules) and some N 2 H 2 (about 100 molecules) are formed during this For comparison, we include a simulation of shock compression of a small sample of the explosive RDX (cyclotrimethylenetrinitramine) at a shock wave speed of 8 km/s. The simulation is performed on a perfect crystal of 5 × 4 × 3 unit cells with 4 molecules in each unit cell. Since no imperfections or air-filled cavities are present in the sample, a lower sensitivity can be expected (compared to a macroscopic inhomogeneous sample). Figure 7 shows (black curve) the temperature as a function of time. The initial shock compression yields a somewhat lower temperature (about 1250 K) than the hydrazine simulation at the same shock speed. However, substantial decomposition after compression results in a steady increase in temperature up to 4000 K after 40 ps.
The reactive molecular dynamics simulations suggest that liquid hydrazine has a lower shock sensitivity than a single crystal of RDX. Furthermore, we distinguish no significant temperature increase during the first 40 ps for a shock speed Vapour pressure as function of temperature for hydrazine [17] of 8 km/s. From these findings, we consider prompt detonation for the 3.5 mm debris particle unlikely. For higher shock pressures, however, simulations indicate a substantial increase in temperature due to chemical reaction approaching that of RDX at 8 km/s.

Hydrazine vapour initiation
The shock wave and detonation properties of hydrazine vapour can be expected to depend on the density of the vapour. We assume that a significant amount of liquid hydrazine is present in the tank at temperatures around and above 300 K. As the tank temperature is increased (e.g., due to radiation exposure from the sun), the vapour pressure (and density) will increase. Figure 9 shows how the vapour pressure of hydrazine depends on the temperature. In this study, we have chosen to more closely examine at three different pressures: 5 bar (443 K), 1 bar (386 K), and 0.2 bar (345 K). For ambient temperature (300 K) and a constant pressure of 5 bar from an inert gas in the tank, we assume the hydrazine to be in liquid form.
Detonation in pure (without oxygen) hydrazine vapour has been reported and investigated in other studies [18,19]. The detonation pressure is many orders of magnitude lower than for a typical condensed explosive, and the shock pressure required for decomposition to proceed is thus much lower. Due to the large mismatch between the shock impedance of the liquid layer (or tank material) and the hydrazine vapour, the shock pressure transmitted into the vapour is much lower than the shock pressure in the condensed material adjacent to the vapour. In this section, our aim is to quantify these phenomena using hydrocode and thermochemical simulations. The critical questions at this point are: Can the pressure wave transmitted to the vapour lead to detonation? and What are the consequences of such a detonation on the tank structure?
In order for a detonation to occur, a sufficiently strong shock wave is required. Furthermore, the wave must have an extension which is comparable to the critical diameter of the reactive material. Pedly et al. [18] report experiments on critical tube diameters as a function of vapour pressure. For example, at an initial pressure of 1 bar, a tube diameter of 23 mm is reported for a sustained detonation to propagate outside the tube (detonation propagates from tube to unconfined environment). At a pressure of 0.2 bar, the critical diameter is increased to about 70 mm. A large debris particle (generating a larger plume) is thus more prone to cause initiation. A smaller fragment may have to travel a longer distance for the expanding material to reach an extension that corresponds to the critical diameter, which will result in a lower shock pressure. Figure 10 shows the pressure field in hydrazine vapour at a pressure of 5 bar at times 0 µs (a), 0.25 µs (b), 5 µs (c), and 10 µs (d) after impact. As the plume is formed, the peak pressure in the shock wave decreases from about 500 bar (0.25 µs) down to 100 bar (10 µs) at a shock wave front width of 35 mm. Simulations are performed for the three different vapour pressures, and the shock wave front is characterised at a position where its width is equal to the critical diameter for each vapour pressure. We conclude that the simulated shock wave pressures are significantly higher (factor of 4-10) than the detonation pressure for the respective vapour pressures. In experimental studies of hydrazine vapour detonation [18,20], hydrazine tubes are initiated by igniting a driver gas (e.g., acetylene-oxygen) which has a detonation pressure similar to that of hydrazine vapour. It is thus concluded that initiation of the hydrazine vapour is likely for all three vapour pressures.

Hydrazine vapour detonation
We use the thermochemical computations to estimate the detonation properties for hydrazine vapour at the chosen pressures. Table 2 summarises the results of these calculations. We note that the calculated detonation velocity and pressure are a bit lower (10-40%) than what has been reported in experimental studies [18,21]. In a next step, we use the thermochemical data as a description of the detonation wave in the finite element code LS-DYNA. This is achieved by adapting the standard model for condensed explosives (referred to as EOS_JWL) to the hydrazine data obtained through thermochemical calcula-tion. Using this approach, a detonation wave with suitable detonation velocity, pressure and particle velocity can be included in a finite element simulation. This model is used to describe the damage process that involves ignition and detonation of the hydrazine vapour after penetration.
The LEO propellant tank is approximated by a spherical Ti alloy structure with a wall thickness of 1 mm. Normal incidence is assumed, and hence, cylindrical symmetry is used with the axis of symmetry coinciding with the path of the particle. Any membrane or inert gas in the tank is neglected. The point of initiation is assumed to be located behind the front side tank wall where the extension of the shock wave has just exceeded the critical diameter (a few cm below entrance hole). Since the shock wave front at the time of initiation is small compared to the tank, a point initiation in a stationary gas is used. Figure 11 shows snapshots from a simulation of hydrazine vapour detonation at 5 bar vapour pressure initiated by the 3.5 mm debris particle. This can be expected to be the worstcase scenario with respect to tank damage. The first image (upper left) shows a spherical detonation wave propagating out from the point of initiation with a detonation pressure of about 90 bar. In the next image (top, middle), the detonation wave hits the tank wall in the area close to the debris entrance hole. The wave reflection results in an increase in pressure to almost 300 bar for a duration of some microseconds. Subsequently, the detonation wave travels down until it hits the bottom of the tank. During the time interval 150-200 µs, the detonation wave-tank wall interaction results in a significantly elevated pressure acting on a large surface area of the tank wall. This pressure (60-150 bar) is of the same order as the burst pressure for this type of tank, and tank rupture can thus be expected from this process alone. In fact, this is also confirmed by the simulation; at about 180 µs after initiation the tank wall is subject to fragmentation. We note that the detonation wave also carries significant momentum which adds a dynamic effect to the wall interaction. These results suggest that for the 5 bar vapour pressure case, initiation followed by tank fragmentation is likely. The last image in Fig. 11 (bottom right) shows how the reflected wave in the combustion gases converges and creates a high-pressure region on the symmetry line of the tank. This converging pressure wave will subsequently be reflected back out and possibly contribute to further fragmentation and acceleration of the tank wall material. This process may be altered by the presence of a membrane (and inert gas) in the tank. Figure 12 shows a snapshot of the pressure field of the detonation wave for the 0.2 bar vapour pressure case. The stagnation pressure in the region of interaction with the wall is 2-5 bar, more than an order of magnitude lower than the burst pressure. We conclude that vapour detonation constitutes a negligible contribution to the damage at 0.2 bar vapour pressure. Fig. 11 Pressure field (Mbar) inside the tank during hydrazine vapour detonation at times 15, 55, 145, 180, and 260 µs after initiation. The debris particle impacts the tank at the top and initiation occurs further down. The initial vapour pressure is 5 bar which can be expected to be the worst case with respect to tank damage. The simulation is subject to rotational symmetry around vertical axis Fig. 12 Snapshot of pressure field 55 µs after initiation for the 0.2 bar vapour pressure case. The thin tank wall is visualised in grey with the entrance hole visible at the top Analysis of the 1 bar vapour pressure simulation shows no sign of wall fragmentation. Plastic deformation in the tank material was clearly observed, suggesting that the 1 bar vapour pressure case corresponds to a transition zone where the impact from the detonation is likely to induce permanent damage. We conclude that in combination with other possible damage mechanisms, such as shock wave and crack propagation in the tank wall (from initial impact) and damage to the rear tank wall due to impacting dispersed material from the plume (e.g., vapour mixed with solid particles), rupture is not unlikely.

Damage assessment
In an experimental study by Shäfer et al. [22], thin-walled cylindrical pressure vessels (containing inert gas) made of aluminium alloy and titanium were impacted by hyperveloc-ity aluminium particles of different sizes. A clear correlation is found between particle kinetic energy and mode of failure. Even if the case presented here has a somewhat different configuration (tank geometry and type of titanium), a comparison is interesting. By calculating the kinetic energy of the 3.5 mm particle (5.9 kJ), we conclude that this projectile scenario would be clearly below the rupture limit reported by Shäfer et al. (20-30 kJ). This suggests that non-reactive damage mechanisms alone would not lead to tank rupture. The 13 mm particle in Sect. 4.1, however, has an energy clearly above the rupture limit.
The present investigation indicates that tank temperature is an important factor with respect to tank damage: an elevated temperature leads to a significant propellant vapour pressure which may cause severe damage upon ignition and detonation. For the 3.5 mm aluminium particle at a speed of 14 km/s relative to the tank, we summarise our main findings as follows: -T > 150 • C: A strong vapour detonation wave is likely to result in violent fragmentation/rupture of the tank. -90 • C < T < 150 • C: This temperature interval can be considered a transition zone where a vapour detonation may cause damage and possible rupture in combination with other damage mechanisms, e.g., debris plume interacting with rear tank wall. -T < 90 • C: The vapour detonation has a negligible effect on the tank structure. The 13 mm debris particle has a kinetic energy of 304 kJ, and rupture is expected irrespectively of propellant vapour pressure. In addition, reactive molecular dynamics simulations suggest that detonation of liquid hydrazine cannot be ruled out. For the 0.5 mm debris particle, no hydrazine reaction or tank fragmentation is predicted.

Conclusions
The present study investigates the effect of hypervelocity impact of spherical aluminium debris particles on a titanium alloy propellant tank. Explicit finite element methods are used to predict the penetration, plume formation, and pressure levels in liquid and gaseous hydrazine behind the tank wall. Reactive molecular dynamics is subsequently used to predict the reactive response in liquid hydrazine at those pressures. For the 3.5 mm debris particle considered in this study, the simulations indicate that an even higher impact pressure is required to cause prompt detonation. Furthermore, we use the finite element method to predict the effect of vapour detonations at different vapour pressures on a simplified tank structure. This assessment indicates that without protective structures, debris particles in the few mm size range at 14 km/s relative velocity may cause fragmentation and rupture when impacting on a hydrazine tank. Furthermore, initiation of hydrazine vapour at high initial pressures may contribute to the damage processes. Thus, shielding from both impacting debris and radiation that may cause heating of the propellant is recommended.