A numerical approach to evaluate temperature-dependent peridynamics damage model for destructive atmospheric entry of spacecraft

Abstract

The evaluation of the on-ground casualty risk assessments due to a controlled or uncontrolled re-entry is highly sensitive to the accurate prediction of fragmentation events during an atmospheric re-entry. The main objective of this study is an investigation into the use of peridynamics (PD) to improve the analysis of fragmentation during atmospheric re-entry with respect to currently adopted semi-empirical approaches. The high temperatures characterising such scenarios may substantially impact fragmentation, which requires appropriate modelling of the damage process within the PD method. The damage models in PD require experimentally determined fracture mechanical properties that are unavailable as a function of temperature. This work proposes a numerical methodology to estimate the PD damage parameters changes with temperature to enable the study of fragmentation during atmospheric re-entry. Initially, tensile-testing simulation experiments are performed in peridynamics to calibrate material parameters for steel and aluminium alloys as a function of temperature. Then, a parametric study is carried out to evaluate the temperature-dependent damage model parameters for the same materials. The applicability of the proposed methodology is showcased using a re-entry test case scenario.

Nomenclature

Abbreviations

0D

Lumped mass thermal model

ATV-1

Jules Verne Automated Transfer Vehicle

D4D

Design for Demise

DKR

Detra-Kemp-Riddel

dof

Degrees of freedom

DLSP

Diagonally Loaded Square Plate

DRAMA

Debris Risk Assessment and Mitigation Analysis

ESA

European Space Agency

FEM

Finite Element Method

FOSTRAD

Free Open Source Tool for Re-entry of Asteroids and Debris

PD

Peridynamics

PMMA

Polymethyl Methacrylate

SAM

Spacecraft Aerothermal Model

SCARAB

Spacecraft Atmospheric Re-entry and Aerothermal Break-up

STL

Stereolithography file

VASP

Vehicle Atmospheric Survivability Project

Symbols

$a,\ c,\ d$

Material constants

$b$

Body force density

$c_p$

Specific heat capacity

$C_P$

Pressure coefficient

$C_{\tau}$

Shear coefficient

$dt$

Timestep

$E$

Young’s modulus

$F$

Forces

$G_0$

Energy release rate

$G$

Shear modulus

$H,\ L,\ R,\ W$

Geometrical parameters

$K$

Stress intensity factor

$m$

Mass

$q,\ Q$

Convective heat flux

$Re$

Reynolds number

$s$

Bond stretch

$St$

Stanton number

$t$

Time

$T$

Temperature

$u$

Displacement

$\alpha$

Energy Accommodation coefficient

$\gamma$

Specific heat ratio

$\delta$

Horizon

$\Delta$

Peridynamic mesh spacing

$\epsilon$

Engineering strain

$\theta$

Flow inclination angle

$\kappa$

Bulk modulus

$\kappa_h$

Hardness modulus

$\mu$

Viscosity

$\nu$

Poisson’s ratio

$\rho$

Density

$\sigma$

Engineering Stress

$\sigma_N$

Normal momentum accommodation coefficient

$\sigma_t$

Tangential momentum accommodation coefficient

$\phi$

Local damage parameter

Subscripts

$\infty$

Free stream properties

$c$

Critical values

$melt$

Properties at melting temperature

$w$

Wall properties

$U$

Properties at ultimate strength

$Y$

Properties at yield strength

1.0 Introduction

Over the years, the increasing space debris population is becoming a high-priority issue to the space-faring communities. European Space Agency (ESA) along with other international space agencies have set up “Debris Mitigation Guidelines“ [1] to restrict the creation of new debris in space by moving satellites to a graveyard orbit at the end of active life or disposing them by re-entering the Earth’s atmosphere. However, atmospheric re-entry may pose a significant risk to the species population and infrastructure on the ground from surviving objects. ESA introduced the Design For Demise (D4D) concepts with its Clean Space initiative to minimise the risk of objects surviving an atmospheric re-entry by promoting complete disintegration during the entry process. The mitigation requirements state that there must be less than 1 in $10^4$ (casualty risk) chance of someone being injured by surviving space debris. Compliance with this requirement can be achieved by performing a controlled de-orbit and re-entry where the surviving debris is made to fall in the open ocean. However, there is an additional increase in mass and cost due to a planned, controlled re-entry prohibiting its usage for several satellite missions. An alternative approach is to ensure a safe uncontrolled re-entry in its end of life disposal strategy. While D4D is a design philosophy and is applicable to all types of re-entries, the current study focuses primarily on the uncontrolled re-entries due to the threat of space debris. Larger spacecraft cannot generally reduce the risk adequately for uncontrolled entries and hence require a controlled re-entry, many medium to small spacecraft have no other choice than to have an uncontrolled re-entry due to financial constraints. Consequently, the latter spacecraft missions may have a casualty risk above $10^{-4}$ which can be reduced by making crucial design changes via different D4D methodologies. A proper re-entry analysis and an estimation of the casualty risk metric [2] is critical to determine if a particular design change reduces the overall risk due to re-entry.

Simulating the re-entry of an object involves several disciplines such as: flight mechanics, trajectory analysis, heat transfer, real gas aerothermochemistry, ablation, hypersonics, subsonics and numerical analysis [3]. Including the high-fidelity aspects of all these phenomena into a re-entry simulation at every point of the trajectory requires a strong coupling between these individual disciplines, and can be computationally challenging and time consuming. Hence, low-fidelity tools have been introduced to mitigate the cost by using analytical approximations of aerodynamic and aerothermodynamic parameters but these tools come at the cost of increased uncertainty. In-general, re-entry analysis codes can be classified into two categories [4], (a) object-oriented [5, 6, 7, 8, 9, 10] and (b) spacecraft-oriented [11, 12]. The object-oriented approach is characterised by approximating the spacecraft model and its components into a set of primitive elements like spheres, cylinders and flat plates, essentially reducing the modelling complexity. These methods usually assume that the primitive elements break off at a fixed altitude (78km) and their trajectories are propagated independently to each other. The individual objects are then deemed to be demised based on material consumption due to melting phenomena. Due to the assumption of a fixed break-up altitude, object-oriented methods generally predict a higher ground risk [4]. On the other hand, spacecraft-oriented tools model the spacecraft to a high granularity level consistent with the real design. The aerodynamic and aerothermodynamic coefficients are calculated for realistic and complex geometries [11] as opposed to the primitive shapes of object-oriented codes. Breakup events are then computed by analysing mechanical loads at pre-defined sections/cuts using simplified breaking stress criteria [13] and accurate thermal load modelling for material ablation.

Novel D4D techniques [14, 15, 16] have demonstrated that the current physical modelling of the spacecraft demise is affected by severe uncertainties in the modelling of the material behaviour or breakup processes making difficult to carefully assess if particular design changes will be effective in reducing the casualty risk. This is mainly due to the fact that it is very difficult to obtain validation data for destructive re-entry experiments. The widely adopted fixed altitude break-up model is based on the Vehicle Atmospheric Survivability Project (VASP) experiments [17, 18] by the United States Air Force office in the early 1970s which mimicked re-entry conditions due to orbital decay and its breakup sequence was analysed. A break-up altitude of around 78km was consistently observed from the experiments. However, only aluminium structures were assumed and this break-up altitude maybe different for different materials. Hence, it becomes essential to understand if the break-up altitude varies because of the variation in materials and the parameters associated with re-entry analysis. It is well known that there are two fragmentation mechanisms, (a) destruction by heating (melting and evaporation) and (b) destruction by forces (fracture between modules). It is generally assumed that the main driver for the break-up process of reentering spacecraft is the thermal load [19]. However, analysis of the MIR space station re-entry revealed that the mechanical break-up of the station modules had taken place at an altitude of 76km, before thermal fragmentation becomes the break-up driver at 69km [20]. As a consequence, there is a need for researching into the dominant fragmentation mechanism of re-entry objects to evaluate and analyse for the correct break-up altitude.

Several experiments have been carried out to analyse and create new fragmentation criteria for re-entry tools [21, 22, 23]. However, these criteria are usually very specific to particular types of connections between different modules of re-entry objects as described in Ref. [24]. Currently, only Spacecraft Aerothermal Model (SAM) re-entry code describes the connection between different modules of the spacecraft [25] and imposes fragmentation criteria to be either purely melt-based, force based, or a temperature-based preset criteria [25]. This study is a first step towards understanding how much of the uncertainty in the structural fragmentation events during an atmospheric re-entry can be reduced by using high-fidelity structural analysis methods.

Traditionally finite element method (FEM)-based tools were used to perform structural analysis during re-entry but an alternative structural analysis tool using peridynamics (PD) methods is more suitable to the current research study. Since the intention is to couple the structural analysis tool with a re-entry code, computational cost and ease of setup become very important parameters for analysis. While FEM tools are known to be well-validated and accurate, they are notoriously difficult to model fracture mechanics, requiring a pre-crack, crack growth laws and surface re-meshing while crack propagation [26]. Also, FEM cannot deal with crack initiation problems thereby requiring a pre-existing crack to be present before the start of a fracture mechanics simulation. PD on the other hand shows promise, alleviating the problems stated above.

PD is a new continuum mechanics formulation developed by Dr.Stewart Silling [27] using non-local formulations. The governing equations are integro-differential in nature and do not contain any spatial derivatives. This enables an inherently easy treatment of discontinuities such as cracks. PD can be a very suitable method to be coupled with a re-entry code owing to the facts that crack initiation and propagation is inherently taken care of. However, PD can be very time consuming to simulate an entire spacecraft because of a few constraints which are discussed in further sections. In the current study, an open source PD software known as “Peridigm”, developed by Sandia National Laboratories [28, 29], is used for all numerical simulations.

The main contribution of this study is three-fold: to introduce the unique capabilities of the PD method to the re-entry community, to compute approximate damage parameters variation as a function of temperatures close to melting, and, to showcase the usability of the PD method using a re-entry test case scenario. In such scenario, high-temperature material phenomena and its fragmentation in the plastic regime requires consideration in the analysis. All materials become weaker when they get hot and their material properties at a particular temperature are dependent on the time of exposure to a particular temperature [30]. The timescale of the flow phenomena during re-entry is in the order of minutes, a typical uncontrolled re-entry trajectory only lasts for around 10–20 minutes [31] from the re-entry interface of 120km. This implies that an object undergoes extreme variation in temperatures in the span of a couple of minutes. To simulate the variation of material properties with respect to temperature, it becomes important to consider the data corresponding to an appropriate exposure time which is suitable to the conditions of re-entry. It would be ideal to consider the material parameter data for the exposure time of 5 minutes or less but the lowest available exposure time in the literature is around 30 minutes. To model crack initiation and propagation using peridynamics, the variation of a few fundamental parameters like critical stress intensity factor $(K_c)$ or critical energy release rate $(G_{\mathrm{oc}})$ (discussed in the next sections) with temperatures close to melting is required. These parameters are obtained experimentally using ASTM standard experiments [32]. Although, the data on $K_c$ of several alloys of materials is abundantly available at room temperature, its subsequent data at higher temperatures is unavailable. To bridge this gap, a methodology in PD is devised as a first assumption. There is also a need to introduce the PD method and its capabilities to the re-entry community as this is the first of such efforts in using PD as a tool to study fragmentation during re-entry. The Jules Verne Automated Transfer Vehicle (ATV-1) [33] was chosen as a reference test case scenario to test the current methodology and the usability of PD in predicting re-entry fragmentation. The structural fragmentation of the solar panels is thought to have occurred before the melt-based fragmentation [12, 33] during ATV-1’s actual re-entry, which is an intriguing aspect of this test case. For a more general use, the solar panel fragmentation is always assumed to occur at 95km altitude for object-oriented simulation purposes. It is interesting to investigate if the PD method can predict the solar panel fragmentation altitude using the current methodology. An in-house, open-source re-entry analysis tool named Free Open Source Tool for Re-entry of Asteroids and Debris (FOSTRAD) [6] is used in the current study. FOSTRAD is a simulation suite that allows for the estimation of aerodynamics and aerothermodynamics of an entry object in a continuum or rarefied hypersonic flow by employing the local panel formulation methods.

The first part of the paper discusses the fundamentals of peridynamics methodology and its equations of motion. The second part of the paper highlights the capabilities of PD method in predicting crack initiation and propagation for static and dynamic loading conditions, and the validation of damage modelling in PD with an experimental study. It is followed by the results and discussion which is further divided into three sections: (1) calibration of materials to the tensile testing data at high temperatures, (2) evaluation of PD damage parameter as a function of temperature for two different materials, and (3) solar panel fragmentation analysis of ATV-1 re-entry. The conclusions are elaborated in the final section.

2.0 Numerical tools

This section presents a brief discussion about the different numerical tools and their formulations. PD theory is used to cover the high-fidelity structural analysis part of the problem, and an in-house re-entry analysis tool named FOSTRAD [6] is used to cover various aspects of atmospheric re-entry.

2.1 Fundamentals of peridynamics

There are two main types of peridynamic formulation as bond-based and state-based PD. Both formulations are set up based on fundamental constraints and conservation laws. The peridynamic domain, being a continuum theory, consists of an infinite amount of material points at x. The material points possess a volume $V_x$ , a mass density $\rho_{x}$ and can be subjected to body force densities b(x, t), displacements u(x, t) or velocities $\dot u(x,t)$ . A brief description of both the above-mentioned formulations is addressed below.

2.1.1 Bond-based PD

The acceleration of any particle at x at time t is given by

(1) \begin{equation} \rho(x) \dot u(x,t) = \int_{H_x}^{} (f(u(x',t) - u(x,t), x'-x))dV_{x'} + b(x,t) \end{equation}

where $H_x$ represents the neighbourhood of x with a horizon of radius $\delta$ , u is the displacement vector, b is the prescribed body force density. The material points within the horizon $(\delta)$ are said to belong to the family $H_x$ and contribute to the net force acting on x (see Fig. 1). The term “bond” refers to the interaction between the material points at x and x ^′. In bond-based PD, the pairwise interactions between material points are equal in magnitude and parallel to the relative position vector in the deformed state (see Fig. 2). Comparing Equation (1) to its classical mechanics counterpart reveals that divergence of the stress tensor is replaced by the integral of the force density [ $(f(u'-u,x',x))$ ] over the neighbourhood $H_x$ . This implies that the function f contains all the constitutive properties of the material. For example, pairwise force density function for a linear elastic isotropic solid is given by

(2) \begin{equation} f(u'-u,x'-x) = cs\frac{y'-y}{|y'-y|} \end{equation}

where $y=x+u$ is the position of the material point in the deformed configuration, c denotes the bond constant, and s is the bond stretch defined by

(3) \begin{equation} s = \frac{|y'-y|-|x'-x|}{|x'-x|} \end{equation}

Figure 1. Peridynamic domain.

Figure 2. Deformation of peridynamic materials points x and x ^′.

The bond constant c for a 3D linear isotropic material can be expressed in terms of material constants of classical continuum mechanics and is given by Equation (4) where, $\kappa$ is the bulk modulus. It is important to note that there are constraints on material parameters for bond-based PD as the peridynamic formulation automatically captures a constant Poisson’s ratio of $0.33$ for 2D and $0.25$ for 3D test cases. Also, plasticity and incompressibility effects cannot be modelled in metals using this formulation. A more detailed explanation for the constraints is provided in Ref. [34].

(4) \begin{equation} c = \frac{18\kappa}{\pi \delta^4} \end{equation}

2.1.2 Ordinary state-based PD

Ordinary state-based PD is a more generalised formulation, which builds upon the assumptions and limitations of bond-based PD. In this formulation the interaction forces do not need to be equal in magnitude and opposite to each other. A force state stores the information about PD forces on a bond and the response of a material point depends on the deformation of all the material points within the horizon $\delta$ . The equation of motion can be written as Equation (5) where, $\underline{T}\langle . \rangle$ represents a force state.

(5) \begin{equation} \rho(x) \dot u(x,t) = \int_{H_x}^{} \left[\underline{T}(x,t)\langle x'-x\rangle - \underline{T}(x',t)\langle x-x'\rangle \right]dV_{x'} + b(x,t) \end{equation}

For a 3D linear elastic isotropic material, the force state can be given by Equation (6),

(6) \begin{equation} \underline{T}(x,t)\langle x'-x \rangle = \left (\frac{2ad\delta}{|x'-x|}\theta(x,t) + cs \right)\frac{y'-y}{|y'-y|} \end{equation}

where a, c and d are material constants given by Equation (7), $\kappa$ is the bulk modulus of elasticity, G is the shear modulus, s is the bond stretch as defined in Equation (3) and $\theta(x,t)$ is the dilatational term, which brings in the effect of neighbouring members (in a $\delta$ size) of the interacting material points on the PD force between them.

(7) \begin{equation} a=\frac{1}{2}\left(\kappa-\frac{5G}{3} \right)\quad c=\frac{15G}{2\pi\delta^5} \quad d=\frac{9}{4\pi\delta^4} \end{equation}

2.1.3 Boundary conditions

The PD equations of motion Equations (1) and (5) are non-linear integro-differential equations and generally require only initial displacement and velocity conditions to obtain a solution. While constraint conditions and external loads are not a necessity, they can be applied to model complex systems. The application of these boundary conditions differ from conventional finite element analysis, as they have to be applied to volumetric material points rather than on surfaces. Figure 3 shows the different regions for the application of boundary conditions. Constraint conditions can be imposed by prescribing displacement and velocity fields in a “fictitious material layer ( $\mathscr{R}_c$ )” along the boundary on a non-zero volume. The extent of this fictitious layer is equal to that of the horizon $\delta$ . External loads do not directly appear in PD equations of motion, and they can be applied as body force density in a “real material later ( $\mathscr{R}_l$ )” along a boundary of non-zero volume. The extent of this layer should be as close to the boundary as possible with depth $\Delta$ . The external loads can be expressed as body force density (b(x,t)) using Equation (8), where F(t) is the external load and $V_\Delta$ is the volume of region $\mathscr{R}_l$ .

(8) \begin{equation} b(x,t) = \frac{1}{V_\Delta}F(t) \end{equation}

Figure 3. Different boundary condition regions of PD domain; R denotes the full domain, $\mathscr{R}_l$ is the region where an external force is applied, B denotes the boundary layer region, and $\mathscr{R}_c$ denotes the region where constraint conditions can be applied.

2.1.4 Damage modelling

As discussed previously, a bond exists between a material point x ^′ with any other material point x in its horizon $\delta$ . The bond breaks when the bond stretch s as given by Equation (3) exceeds a critical value $s_c$ . To describe the connection of a bond, PD uses a history dependent scalar valued function $\mu$ as given by Equation (9). The bond stays intact for $\mu =0$ when the bond stretch is less than the critical stretch, else the bond is broken and the corresponding bond force becomes zero. This implies an increased load on the remaining bonds of a material point which are more likely to break, leading to progressive failure. Cracks naturally appear due to this progressive breaking of bonds and they propagate until complete failure. The local damage at the material point x is defined by the function $\phi$ as given in Equation (10). Function $\phi$ describes the percentage of bonds broken around a material point in its horizon neighbourhood $H_x$ and always lies between 0 and 1. A value of 1 meaning that all the bonds are broken.

(9) \begin{equation} \mu(x'-x,t) = \begin{cases} 1 & if\ \ s < s_c\ \ \forall\ \ t>0 \\ 0 & \mathrm{otherwise} \\ \end{cases} \end{equation}

(10) \begin{equation} \phi(x,t) = 1-\frac{\int_{H_x}^{} \mu(x'-x,t)dV'}{\int_{H_x}^{} dV'} \end{equation}

The variable critical stretch $s_c$ has physical significance as it can be directly related to the critical energy release rate parameter $G_{\mathrm{oc}}$ , an experimentally measurable fracture mechanics quantity. The critical stretch parameter for bond-based PD is given by Equation (11) and for state-based PD is given by Equation (12). These parameters are derived on the basis of Griffith fracture theory, please refer [34] for the complete derivation of Equations (11) and (12). The critical stretch model works well for brittle material and has been validated for several problems [35].

(11) \begin{equation} s_c = \sqrt{\frac{5G_{\mathrm{oc}}}{9\kappa \delta}} \end{equation}

(12) \begin{equation} s_c = \sqrt{\frac{G_{\mathrm{oc}}}{\left[3G + \left(\dfrac{3}{4}\right)^{4}\left(\kappa-\dfrac{5G}{3}\right) \right]\delta}} \end{equation}

where $\kappa$ is the bulk modulus of elasticity and G is the shear modulus of elasticity.

2.1.5 Convergence

Convergence of numerical approximations of the peridynamic equations differs from traditional convergence in the FEM. Several studies have shown that classical continuum mechanics is a subset of peridynamics and that the peridynamic theory converges to the local solution of continuum mechanics for a horizon approaching zero. Zimmermann [36] shows the convergence to the local solution for bond-based PD, Silling and Lehoucq [37] shows that the state-based PD stress tensor reduces to the classical local Piola-Kirchoff stress tensor in the limit of shrinking horizon. Overall, the spacing between the material points, $\Delta$ , and the horizon size, $\delta$ , have significant influence on the numerical results from PD method. It is critical to figure out what these parameters should be to obtain high accuracy while using minimal computational resources.

Multiple specific horizon values are suggested in different publications over the years. In many convergence studies specifically for quasi-static tension [38] and fracture predictions [35, 39] a value of horizon $\delta = 3.015\Delta$ is suggested for homogeneous materials. Hence, in this study, only the variations in mesh spacing ( $\Delta$ ) are considered to comment on the convergence of results.

2.2 Re-entry analysis tool

Free Open Source Tool for Re-entry of Asteroids and Debris (FOSTRAD) is an in-house, medium-fidelity object-oriented code that uses local panel inclination methods to enable the modelling of any arbitrary geometry using small triangular facets (STL files). This enables a rapid computation of the aerodynamics and aerothermodynamic properties of objects during atmospheric re-entry. In addition, FOSTRAD has been modified to incorporate the NRLMSISE-00 atmosphere model [40] for estimating total mass density, atmospheric temperature, and free stream characteristics as a function of altitude for the aerodynamics and aerothermodynamics modules. Furthermore, a $\text{fourth}$ order Runge-Kutta ordinary differential equation solver is used to propagate the trajectory and all object attributes in six degrees of freedom $(6-dof$ ) [41]. This section briefly describes the analytical formulations that are used in FOSTRAD.

2.2.1 Aerodynamics

A surface mesh is obtained by modelling the given object into triangular facets. The pressure and shear stress contribution from each of the facets are computed as a function of local flow inclination angle $(\theta)$ . The resultant force components are computed via the integrals of pressure and shear stress distributions over the surface of the object. A re-entering object experiences different flow regimes, all the way from free-molecular, to transitional and continuum regimes. The aerodynamic contribution of each panel (triangular facet) in the free molecular regime is evaluated using the analytical model of Schaaf and Chambre’s for pressure $(C_p)$ and shear $(C_{\tau})$ coefficients as described in Equations (13) and (14).

(13) \begin{gather} \nonumber C_p = \frac{1}{s^2} \left[ \frac{2-\sigma_N}{\sqrt{\pi}}s\sin{\theta} + \frac{\sigma_N}{2}\sqrt{\frac{T_w}{T_{\infty}}} \right]e^{-(s\sin{\theta})^2} + \\ \nonumber \left[ \frac{(2-\sigma_N)}{s^2}\left((s\sin{\theta})^2+\frac{1}{2}\right) (1+erf(s\sin{\theta})) \right]+\\ \left[ \frac{\sigma_N \sin{\theta}}{2s}\sqrt{\frac{\pi T_w}{T_{\infty}}} (1+erf(s\sin{\theta}))\right]\end{gather}

(14) \begin{gather} C_{\tau} = -\frac{\sigma_t \cos{\theta}}{s\sqrt{\pi}} \left[ e^{-(s\sin{\theta})^2} +\sqrt{\pi}s\sin{\theta}(1+erf(s\sin{\theta})) \right] \end{gather}

(15) \begin{equation} s = \frac{V_\infty}{\sqrt{2RT_\infty}} \end{equation}

where, $\sigma_N$ and $\sigma_t$ are normal and tangential momentum accommodation coefficients. $T_\infty$ is free stream temperature, $T_w$ is the wall temperature, s is the molecular speed ratio given by Equation (15), and $V_\infty$ is the free stream velocity.

The aerodynamic contribution in the continuum regime is computed using modified Newtonian theory as given by Equation (16). $C_p$ is the local pressure coefficient and $C_{p_{\max}}$ is the stagnation point pressure coefficient. The shear stress contribution in the continuum regime is assumed to be zero. The aerodynamic coefficients in the transition regime are calculated with a local-radius-based bridging methodology as defined in Ref. [6, 42].

(16) \begin{equation} C_p = C_{p_{\max}}\sin^{2}{\kern-1pt}(\theta) \end{equation}

2.2.2 Aerothermodynamics

Although different semi-empirical aerothermodynamic models are available in FOSTRAD [6], Detra-Kemp-Riddel (DKR) model is chosen for use in the continuum regime of this study as it is used in SCARAB [13] code, a spacecraft-oriented tool. DKR model uses a Reynolds number formulation to estimate stagnation point aerothermodynamic properties as in Equation (17)

(17) \begin{equation} St = \frac{2.1}{\sqrt{Re_{\infty,0}}} \end{equation}

where St is the stanton number, $Re_{\infty,0}$ and $\mu(T_o)$ are given by Equation (18). $r_N$ is the effective nose radius of the object, $\mu(T_o)$ is stagnation point viscosity and $w = 0.78$ is a temperature exponent of viscosity.

(18) \begin{gather} \nonumber Re_{\infty,0} = \frac{\rho_\infty V_\infty r_N}{\mu(T_o)} \\[5pt] \frac{\mu(T_o)}{\mu(T_\infty)} = \left( \frac{T_o}{T_\infty} \right)^w \end{gather}

The heat transfer in the free molecular regime is given as a function of the local flow inclination angle in Equation (19) where $\alpha$ is the energy accommodation coefficient and $\gamma$ is the specific heat ratio. The aerothermodynamic coefficients in the transition regime are calculated with a local-radius-based bridging methodology as defined in Ref. [6, 42].

(19) \begin{gather} \nonumber q = \frac{\alpha \rho_\infty V_\infty^3}{4s^3\sqrt{\pi}} \bigg[\!\! \left( s^2 +\frac{\gamma}{\gamma-1}-\frac{\gamma+1}{2(\gamma-1)}\left(e^{\frac{T_w}{T_\infty}}\right)\! \right) \\ (e^{-(s\sin{\theta})^2} + \sqrt{\pi}s\sin\theta)(1+erf(s\sin{\theta})) -\frac{1}{2}(e^{-(s\sin{\theta})^2}\bigg] \end{gather}

2.2.3 Thermal analysis

The lumped mass approach (0D) is used to simulate the thermal-ablation of a re-entering object. It is assumed that the temperature of the object is uniformly distributed where multiple objects can have different values of temperature. It is also assumed that the ablation of the object is initiated once the material reaches the melting temperature. The amount of heat remaining after the material reaches the melting temperature $(T_{\mathrm{melt}})$ is given by Equation (20), where $Q_{\mathrm{in}}$ is the total incoming convective heat flux, $Q_{\mathrm{rad}}$ corresponds to the cooling of surface due to radiation, $c_p$ corresponds to the temperature dependent specific heat, dt corresponds to the time step of integration, and, $m_0\ and\ T_0$ correspond to the mass and temperature of the object at previous time-step, respectively. Finally, the mass lost due to ablation is given by Equation (21) where $Q_{\mathrm{melt}}$ is the latent heat of fusion.

(20) \begin{equation} Q_f = m_0c_p \left[ \frac{(Q_{\mathrm{in}}-Q_{\mathrm{rad}})*dt}{m_0c_p} - (T_{\mathrm{melt}}-T_0) \right] \end{equation}

(21) \begin{equation} dm = \frac{-Q_f}{Q_{\mathrm{melt}}} \end{equation}

3.0 Benchmark simulations and validation of the damage modelling

This section provides a few benchmark problems to showcase and highlight the capabilities of PD method which can be useful while simulating the conditions during an atmospheric re-entry. An open source PD software known as “Peridigm”, developed by Sandia National Laboratories [28], is used for all simulations in the current study. Besides the best efforts, the authors did not find any previous literature of using Peridigm [28], code to showcase different novel capabilities of PD theory and hence a few test cases are performed to highlight the same. These test cases are selected from the textbook “Peridynamic theory and its applications” by Erdogan Madenci and Erkan Oterkus [34]. A plate with a circular hole test case is used to showcase the capability of crack initiation and growth under quasi-static loading. A pre-cracked plate test case is used to showcase the capability of crack branching phenomena under dynamic loading conditions. Finally, a diagonally loaded square plate specimen test case is used to validate the damage model used in PD by comparing the results of an experimental study [43]. It is worth noting that comprehensive convergence tests (both $\Delta$ and $\delta$ variations) on these standard test cases have not been performed as the discretisation parameters are taken directly from the references.

3.1 Plate with a circular hole

An isotropic flat plate with a circular hole [34] is subjected to a quasi-static loading and explicit time integration, shown in Fig. 4. This test case is intended to showcase the ability of PD to initiate a crack at the stress concentration zone and its propagation. The geometrical parameters of the test case are as specified in Fig. 4.

(22) \begin{equation} \dot u_y(x,\pm W/2,t) = \pm 2.7541 \times 10^{-7} \end{equation}

Figure 4. Geometry of a plate $(0.050 \times 0.050 \times 0.005\mathrm{m})$ with a circular hole of diameter $D = 10^{-2}\mathrm{m}$ .

The Young’s modulus and density of the material are $E = 192\mathrm{GPa}$ and $\rho = 8,000\mathrm{kg/m}^3$ , respectively. The domain discertisation parameters in PD are shown in Table 1. The top and bottom ends are subjected to a velocity boundary conditions given by Equation (22) over the region $\mathscr{R}_c$ in Fig. 4. Initially, the simulation is run using a quasi-static time integration scheme with a time step of $dt = 1$ s, for a total time of $t_\mathrm{ final} = 1,000$ s in the absence of damage model to obtain the displacement contours due to applied loading. Figure 5 compares the displacement curves along the central X and Y axis at the end of 1,000 time steps, with the ANSYS FEM results [44] and a PD code by Ref. [44]. A close agreement is observed between the current results and the Ref. [34]. Thus, the discretisation parameters shown in Table 1 predict results that are consistent with the literature.

Figure 5. Plate with a circular hole subjected to quasi-static loading when failure is not allowed.

Table 1. PD discretisation parameters

After the above scenario, failure between the material points is allowed by considering a critical stretch value of $s_c = 0.02$ and with boundary condition as $\dot u_y(x,\pm W/2,t) = \pm 2.7541 \times 10^{-3}$ , an explicit solver with a safety factor of $0.9$ is run until a physical time of $t_\mathrm{final} = 0.001$ s. The explicit solver estimates a stable time step based on an unconditional convergence requirement as explained in Ref. [34]. A time step of $\Delta t = 1.399\times 10^{-8}s$ is estimated for the current simulation run. Although there is no pre-defined crack in the Fig. 4, the PD theory predicts the crack initiation at the stress concentration zones. For a plate with a circular hole in uni-axial tension, the stress concentration occurs on the centreline immediately next to the region of hole (see Ref. [45]). Figure 6 clearly shows a self-similar crack growth initiated from the ends of the hole. This is a unique capability of PD theory unlike the other existing FEM techniques which require a pre-defined crack location and a crack-growth law to start the simulation. This particular feature of PD theory can be very useful to study structural integrity during re-entry flows.

Figure 6. Damage contours for plate with a circular hole when failure is allowed.

3.2 Pre-cracked plate

An isotropic flat plate with a pre-defined crack of length $2a = 0.01$ m is shown in Fig. 7. The geometrical, discretisation parameters and material properties are same as that in Section 3.1. The pre-existing crack is defined in the PD simulation by cutting off the interactions between material points across the crack surfaces as described in Ref. [29]. The simulation conditions for the current test case are taken from Ref. [44]. An explicit time integration with a stable time step of $\Delta t = 1.3367\times 10^{-8}$ s using a bond-based PD model is run for a total of 1,300 time steps. A velocity boundary condition of $V_o(t) = 20\mathrm{m/s}$ is applied in the vertical direction onto a region $\mathscr{R}_c$ as shown in Fig. 7. Initially, a simulation is run without considering a damage model to compute the crack opening displacement as shown in Fig. 8. It show cases the cusp like crack-opening behaviour near the crack tip, which is more physically acceptable than the elliptical crack opening prediction of continuum mechanics [27]. Hence, the PD theory captures a meaningful crack opening behaviour.

Figure 7. Geometry of a pre-cracked plate $(0.050 \times 0.050 \times 0.005\mathrm{m})$ with a crack of length $2a = 0.01\mathrm{m}$ .

Figure 8. Crack characteristics.

Now, a simulation is run using the above conditions while considering a critical stretch damage model of $s_c = 0.04472$ . Damage contours are shown at the end of 1,250 time steps in Fig. 9 showing a self-similar crack growth. The crack growth rate can be plotted by considering the position of the crack tip based on the local damage value of any material point exceeding the value of $\phi = 0.5$ , which implies that 50% of the bonds are broken. The crack growth rate obtained from the current study using Peridigm is compared with the curve obtained in Ref. [34] using the PD code in Ref. [44] for the same conditions and the results compare well with each other. From the Fig. 10, crack growth speed can be evaluated to be $1,645\mathrm{m/s}$ from the current study and $1,650\mathrm{m/s}$ from the Ref. [34]. The current crack growth speed is less than that of the Rayleigh wave speed of $2,800\mathrm{m/s}$ , which is an upper limit for the current type of loading [35].

Figure 9. Damage contours at the end of 1,250 time steps under the velocity boundary condition of $\dot{u_o}(t)=20\mathrm{m/s}$ .

Figure 10. Crack growth displacement.

Finally, the applied velocity boundary conditions are increased to a value of $V_o(t) = 70\mathrm{m/s}$ , where the crack growth behaviour changes from self-similar for the earlier conditions to crack branching. The simulation is run for 1,300 time steps just by changing the boundary conditions and keeping the other parameters same as from above. Figure 11 shows the damage contours at the end of 1,000 time steps and a complex crack branching phenomena can be observed. It is worth noting that the crack branching phenomena was achieved without supplying any external branching criteria and thus show casing the advanced capabilities of PD theory in predicting dynamical phenomena.

Figure 11. Damage contours indicating the crack branching phenomena at the end of 1,000 time steps under the velocity boundary condition of $\dot{u_o}(t)=70\mathrm{m/s}$ .

3.3 Diagonally loaded square plate (DLSP)

A diagonally loaded square plate specimen (DLSP) (as shown in Fig. 12) is used to validate the use of critical stretch damage model as a failure criteria. An experimental study using diagonally loaded square plate specimen is carried out by Ayatollahi and Aliha [43] to investigate the effect of mixed mode loading conditions on DLSP specimen. The DLSP specimen has a length $2W = 0.150$ m and a thickness of $t = 0.005$ m with a centre crack of length $2a = 0.045$ m. Two holes of radius $0.004$ m were drilled from a distance of $0.025$ m from the diagonal corners of the DLSP specimen, as shown in Fig. 12. The material is made up of a brittle polymethyl methacrylate (PMMA) material of density $\rho = 1,180\mathrm{kg/m}^3$ , with a modulus of elasticity $E = 2,940\mathrm{MPa}$ and a possion’s ratio of $\nu = 0.38$ . An external tensile load at a constant rate of $1.667\times 10^{-5}\mathrm{m/s}$ is applied on two holes via the loading pins. Also, a critical stress intensity factor for mode-I type of loading on PMMA material is given as $K_{IC} = 1.33\mathrm{MPa}\sqrt{\mathrm{m}}$ . The centre crack propagation paths for various crack angles $\alpha$ are monitored.

Figure 12. Diagonally loaded square plate (DLSP) specimen.

A peridynamic simulation of the above experiment is modelled using bond-based PD to obtain crack propagation paths and compare it with the available experimental data for various crack angles $\alpha$ from Ref. [43]. The domain discretisation parameters are shown in Table 2. A square plate is initially constructed with $150\times 150\times5$ material points of equal spacing and then the material points which lie on the two holes of 4mm radius each are deleted from the resultant domain. A pre-crack of length $2a = 0.045$ m is also modelled in Peridigm. Reference [43] also evaluated the mode-I critical stress intensity factor to be $K_{IC} = 1.33\mathrm{MPa}.\sqrt{\mathrm{m}}$ . According to linear elastic fracture mechanics concepts, the critical energy release rate parameter can be calculated from critical stress intensity factor using the relation Equation (23) where E ^′ is given by Equation (24). These equations can be used to calculate the PD damage model parameter corresponding to the reported $K_{IC}$ value in the experimental work.

(23) \begin{equation} G_{\mathrm{oc}} = \frac{K_\mathrm{IC}^2}{E'} \end{equation}

(24) \begin{equation} E' = \begin{cases} E & \mathrm{plane\ stress}\\[5pt] \dfrac{E}{1-\nu^2} & \mathrm{plane\ strain}\\ \end{cases} \end{equation}

Table 2. PD discretisation parameters

The critical stretch damage model is used to allow for the propagation of cracks in the DLSP specimen and using the Equation (11) gives a value of $s_c = 4.82\times 10^{-3}$ . The loading rate of $1\mathrm{mm/min}$ used in Ref. [43] corresponds to a prescribed displacement boundary condition of $1.6667\times 10^{-5}\mathrm{m/s}$ is applied on a boundary layer of horizon $\delta$ spacing from the hole boundaries. The corresponding simulation is carried out using an explicit time integration with a safety factor of $0.9$ for crack angles $\alpha = 0^0,\ 15^0,\ 30^0,\ $ and $45^0$ . The simulations are carried out until the cracks are fully propagated through out the DLSP specimen. For all fracture orientation angles that are simulated, the crack propagation paths determined from peridynamic simulations and those acquired from experimental results [43] correspond well, as shown in Fig. 13. Crack growth initiation angles [43] are also compared with that of the experimental results. A good comparison is obtained as shown in Fig. 14. This implies that the critical stretch damage model can be used to study crack propagation problems.

Figure 13. Crack path comparison of DLSP experiments [43] and the current peridigm simulations.

Figure 14. Direction of crack growth initiation comparison of DLSP experiments [43] and the current peridigm simulations.

4.0 Results and discussion

4.1 Material calibration

One of the main research contribution of this study is to identify and compute appropriate damage modelling parameters needed for PD to be able to simulate the fragmentation of an object re-entering the Earth’s atmosphere. Materials such as steel and aluminium, which are used in spacecrafts, often have a multitude of alloys that posses drastically different mechanical properties. Before running any peridynamic simulations, there is a need to calibrate the material model to capture the necessary non-linear behaviour of the stress-strain curve.

In this study, an aluminium alloy Al6061-T651 and a 316-annealed stainless steel materials are modelled. The above materials were chosen solely based on the availability of the high-temperature tensile testing data and fracture mechanical parameter data at room temperature. A tensile testing experiment is simulated using Peridigm to calibrate the material model in the following sub-section.

4.1.1 Tensile testing simulation

A tensile testing experiment is modelled and carried out in Peridigm using a standard tensile testing specimen as shown in Fig. 15. The dimensions of the specimen measure $0.1$ m long, $0.012$ m wide, and $0.0025$ m deep as per the ASTM E8M standard. It is important to note that a detailed convergence study (both horizon $\delta$ and material point spacing) on the tensile testing simulation using a state-based PD formulation using Peridigm is already performed in Ref. [46]. The convergence study observed a minimal relative error using five PD material points in the thickness direction and a horizon value that is three times that of the material spacing. Hence, the spacing between the material points in the current simulation is discretised as $0.5\times 10^{-3}\mathrm{m}$ . A computational strain gauge is modelled at the centre of the specimen with dimensions $0.025\mathrm{m}$ long and $0.006\mathrm{m}$ wide. The ends of the specimen are constrained in the x and z directions, and a strain loading rate of $10^{-6}\mathrm{m/s}$ is applied over a horizon distance $\delta$ on the ends. Each of the simulations are run using a quasi-static solver for a physical time of 150s. Since PD theory deals with force densities, the resultant force on the ends of the specimen can be evaluated using Equation (8) and thereby computing the resultant stress on the specimen. The strain in the specimen is directly obtained from the computational strain gauge.

Figure 15. Standard tensile testing specimen as per ASTM E8M standard.

The materials under consideration are calibrated in Peridigm using engineering tensile stress-strain data, for various temperatures from Ref. [47] for $Al6061-T651$ and Ref. [48] for 316-stainless steel. Owing to the current limitations of Peridigm, the non-linear zone of the stress-strain curve and the subsequent plastic zone effects cannot be modelled accurately. So, the two materials are modelled using an isotropic, linear-elastic linear-hardening material model (state-based PD method) with a hardness modulus $\kappa_h$ describing the material response after elastic limit and also requires an input of yield stress $(\sigma_Y)$ . A least-squares approach is used to determine the hardening modulus parameter based on the available experimental tensile testing data. The simulations are carried out for various temperatures by varying the hardness modulus $(\kappa_h)$ and $(\sigma_Y)$ to match the non-linear zone as close to the available data as possible.

The obtained stress-strain curve is plotted against the experimental data of Refs [47] and [48] for Al6061-T651 and 316-stainless steel, respectively. Figures 16 and 17 show the comparison between the Peridigm simulation and the experimental stress-strain data. The calibrated material model compares well with that of the experimental data and the final material model parameters are tabulated in Tables 3 and 4.

Figure 16. Material calibration for Al6061-T651.

Figure 17. Material calibration for 316-annealed steel bar.

Table 3. Elevated temperature material properties of Al6061-T651

Table 4. Elevated temperature material properties of 316-stainless steel

4.2 Damage calibration at high temperatures

To model the structural fragmentation process during an atmospheric re-entry, the fracture mechanical properties of materials at elevated temperatures is necessary. As an initial study, a critical stretch parameter is assumed based on the previous experience and a PD simulation is carried out by applying a preset temperature-dependent loading condition at the boundary to check for failure. If the object does not fail at the preset loading condition, then the critical stretch value is parametrically changed until failure is observed in the geometry under analysis. This approach is currently utilised to calibrate the critical stretch parameter to fail at a particular loading condition provided the material model is calibrated to the tensile testing data. Consequently, this boundary condition should correspond to that of the failure load. But, the data on failure load of a particular material with increasing temperature is not available in the literature and assumptions need to be made for estimating a such a quantity. In this case study, an assumption that the object fails when a load corresponding to that of the ultimate stress $(\sigma_U)$ is applied on the boundary of the geometry is established. Using this methodology the damage parameters at higher temperatures are estimated in the current research.

Under the assumptions of linear elastic fracture mechanics, the parameters describing fracture are material specific and do not depend on the geometrical aspects. Three ideal testing geometries are considered to ensure that the estimated damage parameter is independent of the type of geometry used. Figure 18 shows different geometries that are considered for damage calibration study. Furthermore, PD simulations are run for each specimen at various material point spacings to ensure that the predicted damage parameter is unaffected by discretisation parameters.

Figure 18. Different testing geometries for damage parameter calibration simulations.

Under the assumptions of linear elastic fracture mechanics, the parameters describing fracture are material specific and do not depend on the geometrical aspects. Three ideal testing geometries are considered to ensure that the estimated damage parameter is independent of the type of geometry used. A flat plate of size $(0.1\times 0.1 \times 0.01\text{m})$ , a solid cylinder of size $(0.02\times 0.02\times 0.1 \text{m})$ , and a standard tensile testing specimen (see Fig. 15). Figure 18 shows different geometries that are considered for damage calibration study. Furthermore, PD simulations are run for each specimen at various material point spacings to ensure that the predicted damage parameter is unaffected by discretisation parameters. The ultimate stress variation of the two materials is given in Tables 5 and 6.

Table 5. Ultimate stress variation with temperature of Al6061-T651 [47]

Table 6. Ultimate stress variation with temperature of 316-stainless steel [48]

The simulations have been carried out using the calibrated material model from Section 4.1.1 with an explicit time integration for $15\times 10^{-3}$ s. A stable time step is automatically calculated for each simulation and safety factor of $0.9$ is used on all simulations. The critical stretch damage parameter is parametrically varied until a complete failure is observed due to the uni-axial tensile boundary conditions. Figures 19 and 20 shows the final variation of the critical stretch parameter with temperature for Al6061-T651 alloy and 316-stainless steel alloy, respectively. All three geometries and varying spacing ( $\Delta$ ) parameter give approximately a constant value of critical stretch parameter at a particular temperature. It can also be observed that for a particular geometry, the estimated critical stretch parameter converges as the spacing between material points $(\Delta)$ is reduced. From Figures 19 and 20, a restriction can be imposed on the spacing between the material points to 1mm to further reduce the computational time for additional simulations using this data.

Figure 19. Variation of critical stretch with temperature for Al6061-T651 where fragmentation occurs when ultimate force boundary condition is applied.

Figure 20. Variation of critical stretch with temperature for 316-stainless steel bar where fragmentation occurs when ultimate force boundary condition is applied.

The final variation of PD damage parameter $s_c$ as a function of temperature is shown in Fig. 21 assuming an average of all three geometries at 1mm material point spacing. It is worth noting that these damage parameters are entirely estimated within the PD theory and would need a broader comparison with the experiments performed by the fracture mechanics community. However, the data on experimental parameters like the critical stress intensity factor $(K_\mathrm{IC})$ is only available for room temperature. Table 7 shows the comparison between the critical stretch values obtained from the established values of $K_\mathrm{IC}$ in the literature with the critical stretch values computed from this study. The established values of critical stretch are computed using the relations described in Equations (12) and (23). The error percentage in Table 7 is calculated based on the established $s_c$ values. It is observed that the there is an error of $10.41\%$ for the critical stretch values of the aluminium alloy as opposed to a $59.66\%$ of the steel alloy.

Table 7. Comparison of the critical stretch parameters obtained with the established values of $K_\mathrm{IC}$ from literature

The prediction error observed in Table 7 may be attributed to two aspects, one being the inadequate modelling of material’s plastic regime and the other being the suitability of using critical stretch as a damage criterion for ductile materials. The critical stretch damage model is majorly suitable for brittle materials [35] and may not be entirely suitable to predict fracture in the non-elastic zone. A fracture energy-based failure criterion is necessary to account for the effects of ductile behaviour, which is not available in the current version of Peridigm. Furthermore, the current material model pivots to the yield stress point and the stress-strain response after the elastic zone, is approximated as a linear curve as shown in Figures 16 and 17. The bi-linear material model may directly influence the prediction error. Furthermore, the current data on tensile testing and $K_\mathrm{IC}$ are taken from independent sources. Hence, an in-depth fracture mechanical investigation is necessary to understand various sources of error, which would involve performing tensile and fracture mechanical experiments alongside PD simulations on complementary datasets. This is out of context for the current study, and a future study would need to conduct an experimental campaign to fill in the gaps. However, the current study enables the use of the damage parameter data at high temperatures to model aluminium at a reasonable degree to perform fracture mechanical simulations to study fragmentation during atmospheric re-entry.

Figure 21. Variation of critical stretch with temperature for Al6061-T651 and 316-stainless steel.

4.3 Re-entry test case scenario

As already discussed in the Section 1, the most interesting aspect of the ATV-1 re-entry is the structural fragmentation of the solar panels prior to their complete thermal demise. It is important to note that the actual flight data to mimic a natural re-entry of ATV-1 is unavailable for usage. However, detailed simulations of the ATV-1 re-entry analysis were performed using SCARAB [49] code, a spacecraft-oriented tool, due to a contract for ESA/ESTEC [33, 50]. A high granularity model was used for SCARAB simulations is as seen in Ref. [50]. During its atmospheric re-entry, it has been observed that the joints connecting the solar panels to the main body break-off (all four joints) at around $93.5$ km of altitude and further undergo complete demise due to melt at around 90km of altitude. This is crucial information that will allow us to compare the break-up altitude of these joints with the current methodology.

4.3.1 Trajectory simulation

The ATV-1 test case is modelled as nine objects involving the main body, four solar panels and four joints connecting the solar panels to the main body, as shown in Fig. 22. The test case is modelled using primitive objects because of the unavailability of detailed geometrical specifications of various components. The main body of the ATV-1 is approximated to that of a hollow cylinder with thickness corresponding to a mass at re-entry interface (33,50). The pressurised shell of the actual ATV-1 and the internal structure is made up of Al2219 and $Al6061-T6$ material, respectively. However, the current model is assumed to be entirely made of homogeneous $Al6061-T651$ material to utilise the simulated results in PD. The geometrical specifications used in the current study are elaborated in Table 8. It is worth noting that there is no information about the mass of the joints and therefore its thickness is assumed to be $0.01$ m, which corresponds to a mass of $25.53$ kg, for the initial set of re-entry simulations.

Figure 22. Object-oriented model of the Jules Verne Automated Transfer Vehicle (ATV-1).

Table 8. Geometrical parameters of the ATV-1

The initial conditions at the 120km re-entry altitude in the geodetic coordinate reference frame are derived from the references [33, 51] and are tabulated in Table 9. The re-entry trajectory of the ATV-1 is simulated using the in-house developed FOSTRAD code, ESA’s DRAMA tool [7] and compared with that obtained from the results using SCARAB [51]. Figure 23 shows the comparison between the main body trajectories of the ATV-1 re-entry simulations. Firstly, the trajectories from FOSTRAD and DRAMA are in strong agreement with each other. This is understandable as both the tools use the same geometrical model as in Fig. 22. The re-entry trajectory from the current simulation has a good agreement with that of the SCARAB simulation until around 80km altitude after which the trajectory diverges from the reference. This may be because of high-granularity modelling and high-fidelity nature of SCARAB simulation. Nonetheless, the trajectory comparison is deemed satisfactory for the current study as we are only interested in estimating the structural fragmentation altitude of the solar panels. This usually occurs before reaching 80km of altitude.

Figure 23. Comparison of the re-entry trajectory of ATV-1 from FOSTRAD, DRAMA and SCARAB simulations.

Table 9. Trajectory initial conditions at the re-entry interface of 120km altitude

4.3.2 Determination of the fragmentation altitude

The solar panel fragmentation altitude during the ATV-1 re-entry is determined by performing a preliminary coupling simulation using FOSTRAD and Peridigm tools. The coupled simulation starts by specifying the geometry files of the involved components, the material and geometrical characteristics such as the thickness of the individual components and initial conditions in the geodetic coordinate reference frame, at the 120km re-entry interface Table 9. The in-house re-entry analysis code FOSTRAD starts solving the trajectory and monitors the integrated aerothermodynamic loads on the flagged objects (main body and solar panels). The FOSTRAD tool continues its trajectory analysis until the loads on the objects reach a heuristically determined threshold value. The threshold value is determined based on the prior experience with the fracture mechanical PD simulation on the required joints. This assumption reduces the overall computational resources and is used as a switching criterion to initiate PD simulations using Peridigm.

After this point at every time iteration, a peridynamics simulation is called forward on the joints connecting the flagged objects. The resultant aerodynamic forces and moments are transferred to the ends of the joint as boundary conditions and the material properties corresponding to the surface temperature, from the lumped mass approach, are selected. A PD simulation is then performed with an explicit time integration using a time step of $10^{-8}$ s for a total time of the order of O(1)s and the joint is discretised using a material spacing of $10^{-3}$ m. This timescale in PD simulation is chosen to be on par with the time step for the trajectory integration. Then, the damage contours from the PD simulation result are analysed for crack characteristics. A crack is said to have formed if the damage value, using Equation (10), on the material point exceeds $0.5$ . The damage contours are then plotted to determine the length of the crack in comparison with the characteristic length to determine if the joint fragmentation has occurred. The joint is assumed to fail once the crack length is same as the diameter of the joint. If the joint does not fragment then the FOSTRAD simulation is carried on to the next time step in the trajectory solver and the process continues. If the joint fragmentation has occurred then the break-up altitude is noted and the simulation continues till the altitude reaches 0 or until all the objects are completely demised due to melting.

A re-entry analysis of the ATV-1 is performed using the above methodology. From Fig. 24, it can be observed that the temperature of the solar panels reach close to the melting point of the material. Note that the material model and damage parameters are only calibrated for maximum temperature of $500^{\circ}\mathrm{C}$ for $Al6061-T651$ in the previous sections, and the melting temperature for the same is $680^{\circ}\mathrm{C}$ . In the current test case scenario, material and damage parameters are assumed to be corresponding to that at $500^{\circ}\mathrm{C}$ for all temperatures above this value. Also, for this preliminary study purposes, it is assumed that the joint connecting the solar panel to the main body retains the temperature of the solar panel while performing PD simulations. The loading history on the centre of mass of the solar panels with respect to the altitude are shown in Fig. 25. The oscillations in the loading history can be observed from Fig. 25 as a consequence of a very high initial pitch axis spin rate. As discussed previously, a heuristic resultant force value of 600N is chosen as a switching criterion in the current study to enable the start of PD simulations.

Figure 24. Temperature history of the solar panels from FOSTRAD simulation.

Figure 25. Loading history on the solar panels from FOSTRAD simulation.

From the results in Section 4.2, the joint connecting the solar panel to the main body is discretised with a material spacing on $10^{-3}$ m. The loads from the centre of mass of the solar panel in Fig. 25 are translated to one end of the joint and a corresponding moment is also applied as the boundary condition. The other end of the joint is constrained to be fixed as it is connected to the main body which has higher inertia. As this is a preliminary study into the usage of PD approach for re-entry applications, only aerodynamic loads are provided as boundary conditions, and inertial forces are completely neglected in the current research work. A PD simulation is then performed with an explicit time integration using a time step of $10^{-8}$ s for a total time of the order of O(1)s. Figure 26 shows the damage contours on the joint at an altitude of $89.66$ km. It can easily be seen that the length of the crack is comparable to the diameter of the joint, indicating that the joint has failed. From the current methodology, all the four joints connecting the solar panels are found to demise at $89.66$ km altitude as opposed to a $93.5$ km altitude prediction from SCARAB simulations [33]. The detached solar panels further undergo complete demise at $86.2$ km altitude as opposed to a 90km altitude from SCARAB simulations.

Although the current methodology predicts break-up and thermal demise altitudes close to the reference data, it is crucial to highlight the sources of error that may influence the results in reality. These errors could be due to the limitations in the geometrical and physical model of the ATV-1 (see Section 4.3.1), and the errors in the damage parameter calibration at higher temperatures (see Section 4.2). It is worth noting that the material and damage model parameters for $Al6061-T651$ could only be calibrated until a maximum temperature of $500^{\circ}\mathrm{C}$ due to the lack of availability of experimental data. From Fig. 24, it can be observed that the temperature of the solar panels exceeds well above $500^{\circ}\mathrm{C}$ before the fragmentation event. Since the material gets weaker as the surface temperature increases, a higher break-up altitude can be envisaged in reality. Also, it is important to note that the thickness of the joints appears to be important in determining the altitude of solar panel fragmentation. As a result, a parametric study is carried out using the current methodology for various thicknesses of joints linking solar panels to the main body, as shown in Table 10. It can be seen that all four joints connecting the solar panels to the main body break up at the same altitude for a given joint thickness. The break-up altitude reduces as the joint thickness decreases, from $89.66$ km for the highest value of thickness to $91.19$ km for the lowest value of thickness. More importantly, the predictions using the current methodology are not far off from the actual SCARAB predictions. This provides some reasonable level of confidence in using the PD methodology as a tool to understand the fragmentation process during atmospheric re-entry.

Table 10. Break-up altitude of solar panels as a function of thickness of the joint

Figure 26. Peridynamics simulation showing the crack formation on one of the joint connecting the solar panel to the main body.

5.0 Conclusion

Peridynamic (PD) theory is used to determine the damage model parameters needed to mimic the fragmentation of re-entering objects subjected to substantial temperature changes. To demonstrate the usability of PD methodology, a realistic re-entry test case scenario of the Jules Verne Automated Transfer Vehicle (ATV-1) is simulated using an in-house re-entry analysis code named FOSTRAD. The open-source PD software known as Peridigm is used to develop the proposed numerical approach. The requirement for the fracture mechanical data as a function of temperature is highlighted and a methodology based on the ultimate stress criteria is established to evaluate such parameters.

Material calibration was performed for Al6061-T651 and 316-annealed stainless steel by modelling a tensile testing experiment and the results are then compared with the experimental data. Three different geometries were considered to perform a parametric study to evaluate the peridynamic damage modelling parameters as a function of temperature. The results of the peridynamic damage modelling parameters were then compared with the established literature data at room temperature for the two materials under study. It was observed that the results obtained using the current methodology provide satisfactory results for the material made of aluminium while discrepancies were observed in the case of steel, at a room temperature. This may be attributed to the inadequate modelling of the material’s plastic regime with a linear approximation after the elastic limit. While, more advanced material models considering a piece-wise linear modelling of the plastic regime are available for a general peridynamic implementation, they are not available in the open-source version of Peridigm. The results presented here are believed to be the best that can be achieved given the current limitations of the open-source tools. Improvements in the current predictions can be achieved with the use of more advanced material models capturing ductile behaviour by modelling the effects of plasticity. These advanced models are already conceptually proven and demonstrated on the non-open source version of Peridigm. It is envisaged that a parallel experimental and numerical campaign would be required to fill the gaps between the need for these parameters and the ability of Peridigm to simulate certain physics essential for high-temperature fragmentation analysis.

A re-entry test case scenario of the ATV-1 is used to test the current methodology and the usability of PD in predicting the solar panel fragmentation altitude. ATV-1 model is approximated using primitive geometries due to the unavailability of detailed geometrical specifications. The re-entry trajectory from the current simulation is in good agreement with the data in literature until 80km of altitude, after which the trajectory diverges from the reference data. Nonetheless, this is satisfactory as the current interest is in determining solar panel fragmentation altitude, which occurs before 80km. During the re-entry trajectory, the aerodynamic and thermal loads on the solar panels were continually monitored, and a heuristically derived threshold value was chosen as a switching condition to allow Peridigm simulations on the joint linking the solar panel to the main body. The resultant aerodynamic forces and moments are transferred to the ends of the joint as boundary conditions, and the calibrated material and damage model is used to perform Peridigm simulations. All the four joints connecting the solar panels to the main body are found to demise at $89.66$ km as opposed to a $93.5$ km altitude prediction from the reference data. The detached solar panels further undergo a complete thermal demise at $86.2$ km as opposed to a 90km prediction from the reference data. These errors can be attributed to constraints in the geometrical and physical models of ATV-1 and material modelling limitations at higher temperatures. Nonetheless, the predictions from the current methodology are not far off from the reference data. The results obtained in the current study constitute a first critical step for the re-entry community to adopt peridynamics as an approach to perform fracture mechanical simulations at temperatures close to melting.