High fidelity model for the atmospheric re-entry of CubeSats equipped with the Drag De-Orbit Device

The use of CubeSat-like small satellites is growing exponentially nowadays, pushing towards missions of increased complexity, including Earth imaging, commercial communications and astronomical observations. As such, they might require components that may survive the re-entry conditions and reach the ground, posing risks for population and properties, or that are intended to be retrieved. The possibility of demise and ground impact poses many challenges from the modeling standpoint because of the uncertainties associated with both the aero-and the aerothermo-dynamic models of the spacecraft. Several formulations and correlations can be found in the literature. Most of them are derived in dated and difficult-to-access papers and technical reports. This paper collects all the necessary and sufficient models, laws and data to describe in a comprehensive way the re-entry of small satellites. They are presented in an organized fashion, with uniform nomenclature and consistent assumptions in order to provide the smallsats scientific community with a smallsats specific, easy-to-understand and rapid-to-implement tool. Furthermore, the paper originally presents an approximated aero-and aerothermo-dynamic model of the Drag De-Orbit Device, a recently developed drag modulation device for drag-based controlled re-entry of large CubeSats.


Introduction
In recent years, the interest and demand for small satellites have grown exponentially.CubeSat-like [1] nanosatellites are frequently chosen and built by Universities and small non-accademic organizations to exploit their characteristics of flexibility and affordability for onorbit technology demonstration and experimentation [2].The end-oflife design for this type of spacecraft is often approximate or entirely neglected for two main reasons: 1) the materials used to build the spacecraft and its payload are usually designed for demise during the re-entry [1] and 2) a comprehensive and generally adopted re-entry safety policy for CubeSat missions does not exist [3,4].Controlling the decay of a spacecraft through drag modulation has been proven in the past [5][6][7] to be a very effective way to remove the spacecraft from crowded Low Earth Orbits (LEO) after the completion of their mission.This, in turn, provides several benefits from the operational and mission objectives standpoints: 1) it has beneficial consequences on safety of population through the active control of the deorbit point [3], 2) it assures compliance with the 25 years orbital lifetime debris mitigation constraint [4], 3) it provides a certain capability of spacecraft manoeuvering, potentially useful also for on orbit collision avoidance purposes [8], and 4) it allows retrieval of parts/data designed to survive the re-entry.
In this paper we are specifically interested in the Drag De-Orbit Device (D3) [9], shown in Fig. 1.It has been conceptualized in-house and designed as a stand-alone drag device capable of on board guidance computation and tracking.This is achieved through the modulation of the spacecraft's ballistic coefficient to autonomously maneuver the CubeSat to the desired de-orbit location and to provide 3-axis stabilization.The spacecraft's ballistic coefficient is modified deploying or retracting four measuring-tape-style booms, each of which is 4 cm wide when flat, 3.7 m long when fully deployed, and inclined at 20°relative to the rear face of the satellite.
Regardless of the specific strategy, any controlled deorbitation design study, including the D3 controlled decay, has to inevitably deal with the simulation of the atmospheric re-entry trajectory, including the spacecraft fragmentation.Some software have been developed in the past by space agencies and research institutes specifically to address the re-entry fragmentation of space vehicles.Well-known are the NASA's DAS (Debris Assessment Software) [10] and ORSAT (Object Reentry Survival Analysis Tool) [11,12] and ESA's SCARAB (Spacecraft Atmospheric Re-entry and Aerothermal Breakup) [13]  Analysis) [14].Lips and Fritsche [15] give an exhaustive overview.
To characterize the evolution of the spacecraft dynamics from orbital altitude to the ground impact or to a potential fragmentation in the atmosphere, all these software include three main building blocks: the aerodynamic model, the aerothermodynamic model and the dynamic model of the spacecraft [15].More specifically, the aerodynamic model gives an estimation of the drag force acting on the vehicle, characterized as the exchange of momentum between it and the surrounding flow [16].The aerothermodynamic model forecasts the heat power that enters into the structure due to the fact that air flows at high speed about the spacecraft.The friction between the fluid filaments as well as the compression in the neighbourhood of the stagnation point convert the kinetic energy of motion into heat that warms up the flow and consequently the spacecraft [17].
The aerodynamic model and the aerothermodynamic model are nested inside the dynamic model and ultimately provide it with the estimation of two key parameters: the drag coefficient and the heat power at the wall, respectively.Finally, the dynamic model is in charge to provide the trajectory followed by the spacecraft as position and velocity in time through the integration of the equations of motion [18].
Being able to perform an accurate aero-and aerothermodynamic analysis of a spacecraft re-entering the atmosphere is extremely important either to design its thermal shield, in case it has to withstand the re-entry conditions, or to estimate the debris field survivability and dispersion, in case of a destructive re-entry.This type of analysis is however rather complicated by the different atmospheric regions and velocity regimes in which the spacecraft travels and the uncertainty that still exists on several phenomena of dynamic and thermodynamic interaction between the spacecraft structure and the flow field [16,19,20].
Semi-analytical laws and correlations, developed in the literature [21][22][23][24][25], can be used to estimate the key parameters of the two models in order to avoid very time-consuming CFD-based calculations.Nevertheless, the choice of the specific mathematical law is not straightforward.Firstly, because there exists a myriad of available laws and correlations.See for example [17,20,23,26,27].Most of them are derived or described in dated and difficult-to-access papers and technical reports, each using a completely different nomenclature and notation.Furthermore, they are based on specific assumptions, which are often not mutually compatible either one another or with the problem of interest.Secondly, because the choice depends on the flowfield regimes, distinguished through both the Knudsen and the Mach numbers, as well as on the geometric characteristics of the stagnation point, which have to be coherently matched with the assumptions of the selected law.
The first contribution of this paper to the State of the Art is the collection of all the necessary and sufficient models, laws and data to describe in a comprehensive way the re-entry of small satellites.The selected laws are in accordance in terms of assumptions and degree of approximation.They are presented in an organized fashion and with an uniform nomenclature in order to provide the smallsats scientific community with an easy-to-understand and rapid-to-implement tool.Table 1 gives a summary of all the selected correlations.The paper includes also the thermodynamic characterization of the flowfield and guides step by step the reader in the computation of all the required variables.
As a further contribution, the paper originally presents an approximated model, based on existing correlations collected in Table 3, of the D3 system to account for the heat that enters into the deployed booms and to estimate the corresponding bending and demise altitudes.
Finally, the work presented in this paper can be regarded as first building block for an accurate re-entry safety analysis of particular smallsats missions.We are referring here to controlled re-entry scenarios having two peculiarities: 1) targeting smaller-than-usual flight path angles at the deorbit point or, in other words, performing very shallow re-entries, and 2) equipped with components that might withstand the re-entry and reach the ground.The first one is typically the case of CubeSats equipped with the D3 system.Indeed, CubeSats usually are enough small to not have available or necessitate an impulsive propulsive system, so they choose to exploit the D3 system to target the deorbit point.It progressively reduces the semi-major axis accordingly to the available control effort.At the deorbit point (about 120 km altitude), the orbit is still nearly circular, which implies approximately zero flight path angle.This type of re-entry trajectory spends a much longer time inside the atmosphere if compared with the classic steep re-entry.Consequently, the fragments resulting from the spacecraft erosion and melting process spread out extensively inside the atmosphere.If these fragments do not demise, very large areas on ground may be interested by their falling, posing thus risks for population and properties.For this type of missions, a robust and reliable reentry safety analysis is required.It is the topic which reference [28] is focused on and we suggest an interested reader to refer to it for further details.The high fidelity model of re-entry presented in this paper is the first fundamental step in Ref. [28] to statistically characterize the input space of the associated uncertainties.Legge's bridging law [23] Continuum regime < Kn 0.01 Hypersonic regime > Ma 10 Blunt nose Modified Newton Law [22] Fay and Riddel's correl [24].Sharp nose Newton Law [16,32] Eckert's model [25] Hypersonic-supersonic transition Blunt nose Sigmoid bridging formula Fay and Riddel's correl [24].if > Ma 6 or else no heat [13] Sharp nose Sigmoid bridging formula Eckert's model [25] if > Ma 6 or else no heat [13] Low Mach number < Ma 2 Hulburt's data [31] No heat [13] This paper is organized as follows.Section 2 gives a first distinction of the flowfield regimes accordingly to the Knudsen number.Both the aerodynamic and the aerothermodynamic models will be treated accordingly to this first distinction in sections 3 and 4, respectively.Subsequently, the dynamic model is presented in section 5, including a simplified approach to account for the melting of the spacecraft structure.Section 7 focuses on the description of the model used to characterize the booms of the D3 system along the re-entry.A series of test cases simulations with the sole CubeSat in different attitude modes and geometries is given in section 6.In this section, a results comparison is also provided of some re-entry trajectories of box shaped objects presented by Lips et al. [29], run with ORSAT and SCARAB software.Finally, the same study cases in sec.6 are recalled in section 8, but considering the CubeSat equipped with the D3 system.

Flow regimes
Let us define the Knudsen number as: where MFP is the mean free path between successive collisions of air molecules and L c is the flowfield characteristic length, approximated with the body characteristic length.The mean free path can be expressed as [19]: is the effective diameter of the gas particles and n d is the number density that can be expressed as ratio between the air density and the molecule mass.The latter is considered constant and equal to kg kmol 28.9 / . The Knudsen number can be used to distinguish among three flowfield regimes and, as suggested by Wilmoth et al. [30], the following boundaries are considered: • Free molecular regime: > Kn 10 • Transition regime: < < Kn 0.01 10 The Free molecular regime is characterized by a large mean free path and consequently by a long relaxation time.The continuum hypothesis does not hold and it is necessary to consider the state of the single particles and their interaction with other particles and with the boundaries.On the contrary, in the Continuum regime the flow has a very short relaxation time and the macroscopic properties can be considered to vary continuously.Between these two extremes in the Transition regime, both the molecule-surface collisions and the intermolecular forces are important.In the following sections, the aerodynamic and aerothermodynamic models will be treated accordingly to this distinction in the flowfield regime.

Free molecular regime
Assuming the Maxwellian velocity distribution of the particles, the Schaaf and Chambre's analytical model [21] can be used to get a compromise between the perfectly specular and perfectly diffuse reflection model.Breaking the spacecraft's structure down into a series of panels exposed to the incoming flow, the pressure and shear stress coefficients, c Pi and c i respectively, on i-th surface of the body can be expressed as function of the local inclination angle i between the surface panel and the freestream as: (3) where N and T are the normal and the tangential momentum accommodation coefficients, dependent on the material of the body and on the type of gas.A good agreement with experimental data is reported in Hulburt [31], considering 1

N
for interaction of air with most metallic surfaces.Also T is expected to be very close to unity.Here they are considered both equal to 0.9 for Aluminum with air.s is the freestream molecular speed ratio: where = R J kgK 287 /( ) is the ideal gas constant for air.T is the tem- perature of freestream air and T w is the wall temperature of the external structure of the spacecraft.T w is assumed to be available along the numerical integration of the dynamic model (see section 5), being an independent variable of the problem.Obviously, the above relations (Eqs.( 3) and ( 4)) are valid only for panels exposed to the flow, that is Pi is forced to be zero for 0 i and c i for < 0.Thus, the drag coef- ficient for the Free molecular regime can be computed summing over all the body panels:  Deployed Modified Newton Law [22] Fay and Riddel's correl [24].if > Ma 6 or else no heat [13]

Bent
Li and Nagamatsu theory [47] Fay and Riddel's correl [24].if > Ma 6 or else no heat [13] Low Mach number Not of interest No heat [13] where V ˆis the freestream velocity unit vector expressed in the body reference frame, n ˆi is the panel outward unit vector and i is the unit vector tangent to the panel within the plane individuated by V ˆand n ˆi.The inclination angle i between the panel and the freestream is computed as: To model a CubeSat, six panels are considered corresponding to the six facets of the rectangular prism.So each surface dimension S i is computed accordingly to the type of CubeSat.Eight different types of CubeSats are analyzed and a three digit numerical code is used to categorize them indicating the length of each of the three dimensions, length L, width W and height H, respectively, in 1U units (approximately 10 cm).Table 2 lists the dimensions of all the analyzed CubeSats configurations.The reference area in Eq. ( 7) S Ref is always computed as the frontal area ( × H W ) of the CubeSat.

Transition regime
The Transition regime is still nowadays not well understood as the Free molecular and Continuum regimes because of the complex motion of the particles which start to be affected by the intermolecular forces.Usually, global transitional flow bridging relations are used.Here, the bridging formula proposed by Wilmoth et al. [30] is suggested:

Continuum regime
Within the Continuum regime, an additional separation of flowfield regimes must be done accordingly to the Mach number Ma: • Fully established hypersonic regime: > Ma 10; • Hypersonic/supersonic transition regime: < < Ma 2 10; • Supersonic, transonic and subsonic regime: Accordingly to this subdivision, the specific correlation is selected.

Fully established hypersonic regime
Despite the inherently nonlinear behaviour of hypersonic flows, a local surface inclination method can be used to estimate the pressure distribution on hypersonic bodies [16].Indeed, the famous "sinesquared law" postulated by Isaac Newton in his Principia in 1687 [32], even though it has been proven to not be accurate for low-speed flow, gives an excellent prediction of the pressure coefficient at hypersonic regime.It entirely neglects the random motion of the particles associated to the static pressure and so they are considered to move in a perfectly rectilinear stream without any interaction with each other.Thus, the pressure exerted on the surface originates only from the exchange of normal momentum with the surface with the total conservation of the tangential momentum.Because of the extremely high speed of hypersonic flow, the reality closely approaches this ideal situation.The Newton's Law reads: where i is the same inclination angle introduced above.Lees proposes [22] a Modified Newton Law (MNL) very close to Newton's original one but proven to be more accurate for blunt body flowfield.It reads: where c PMax is the maximum value of the pressure coefficient evaluated at the stagnation point: The total pressure at the stagnation point P t2 is computed as it will be explained later in sec.4.3.In this work, both laws are used, accordingly to the attitude of the spacecraft during the re-entry.To simplify the problem, four different attitude modes are considered, as sketched in Fig. 2.
When the spacecraft is in face pointing mode, the nose is very blunt, nearly flat, so the MNL is more appropriate.On the contrary, for edge and corner pointing modes we consider a rather sharp nose and use the original Newton Law.Finally, the drag coefficient is obtained summing over all the panels as reported in eq. ( 6) but neglecting the shear contribution as prescribed by the Newton theory.
To compute the drag coefficient in tumbling mode, a very simple formula, averaging among the other attitude modes, is proposed.Using an averaging approach is a common practise in literature [12,33].The main advantage of this approach is that it does not require to know or compute the attitude of the spacecraft.Thus, as described in sec.5, the spacecraft can be approximated as a point mass, having only three degrees of freedom.1  The simulation results are thus strongly simplified and speeded-up because the attitude propagator is not necessary.The main assumption introduced by this averaging approach is that all the different attitudes modes are equally likely to occur and the spacecraft spends the same amount of time in each mode.As a consequence, its main limitation is that it is not able to foresee potential spinning stabilization and it neglects a possible rotationtranslation coupling.However, some authors [34] found a good accuracy of these simplified averaging approaches and, since their speed and easy implementation, they are considered here as the most suitable for this work.
Generally, the CubeSat may have the three dimensions, width W, height H and length L, all different.In this case, the drag coefficient changes as function of which specific face or edge of the CubeSat is actually pointing to the flowfield.Note that, in a rectangular prism, there are always 2 instances of each distinct face, 4 instances of each distinct edge and the 8 corners are all the same.So we distinguish for face pointing among three distinct coefficients: .Therefore, the drag coefficient for the tumbling mode is computed with a simple average among the several possible attitude modes but weighted on the number of faces, edges and corners, as: where 26 is the total number of faces, edges and corners.
Hypersonic/supersonic transition regime.This transition flowfield regime encompasses a wide Mach range because the characteristic phenomena of the hypersonic flowfield (thin shock layer, viscous interaction, high temperature flows, …) become progressively less important as the Mach number decreases and there is not a precise demarcation limit [16].To mathematically average in this Mach range the drag coefficient with a continuous and differentiable function, the Sigmoid bridging formula: is proposed, where the parameter t is computed through a linear interpolation as: Sup Hyp Sup (17) where = Ma 2 Sup and = Ma 10 Hyp are the selected boundaries of fully established supersonic and hypersonic regime.Thus, the drag coefficient is: where c D Sup and c D Hyp are the drag coefficients computed at the local Mach number Ma but considering as if the flowfield was fully established supersonic or hypersonic, respectively.Supersonic, transonic and subsonic regime.These relatively low-speed flowfield regimes are of side interest for the purpose of this work since it is very unlikely that the spacecraft reaches these regimes still intact.They occur at the very end of the re-entry trajectory, which is nearly vertical, the altitude is almost approaching zero, generally lower than 20 km, and the velocity is the terminal velocity (approximately 50-60 m/s).A simple fitting of experimental data is proposed to directly estimate the drag coefficient, taken from Hoerner [35] and reported in Fig. 3.

Aerothermodynamic model
Three different sources/sinks of heat coming into the structure are considered in the model: heating by convection with the air flow, heating by radiation with the air flow (only in Continuum regime) and cooling by radiation with the environment. 2The first two contributions are collected in the variable Q ExchAir representing the total thermal ex- change with the air flow.The paragraphs below are dedicated to describe in detail how to estimate this term.The third contribution, due to the radiation to the environment and valid for any flowfield regime, is computed as: where = 0.18

Al
is the emissivity dependent on the spacecraft body material, Aluminum in our case, = × W m K 5.67 10 /( ) is the Stefan-Boltzmann constant and T w is, as above, the wall temperature of the external structure of the spacecraft.S Ext is the external surface of the spacecraft, calculated as: where W, H, L are width, height and length of the CubeSat.So, the total heat power entering into the spacecraft structure is:

Free molecular regime
As for the aerodynamic model, the Schaaf and Chambre's analytical model [21] is able to provide also an estimation of the thermal power per unit area reaching the wall.The analytical formulation as function of the local inclination angle i is: where γ is the specific heat ratio considered constant and equal to 1.4 and a c is the thermal accommodation coefficient assumed equal to 1 for metallic materials and air.Some data can be found in Hulburt [31].
Then, the total thermal power entering into the structure due to the thermal exchange with the air flow is simply computed as a discrete summation on the facets of the CubeSat:

Transition regime
The bridging formula used in this work to estimate the aerothermodynamics in the Transition regime follows the approach of Legge [23].Accordingly:

Continuum regime
Within the Continuum regime, the hypersonic regime is extremely important in the overall heat budget that enters inside the body because it contributes for almost half of it.It is therefore treated in detail considering the hypothesis of chemically reacting gas mixture in thermodynamic and chemical equilibrium, after the bow shock wave that develops in front of the spacecraft, and within the boundary layer.Since the extreme viscous dissipation that occurs within the boundary layer and the sudden compression across the shock wave, a sharp increase of the flow temperature is expected, reaching values sufficient to cause dissociation and ionization of the gas.Consequently, the hypothesis of calorically perfect gas is not admissible for this work and so two iterative processes are set up: the first one to compute the air properties across the shock wave and the second one to get the gas mixture composition.In the following paragraphs, it will be explained how to thermodynamically characterize the flow, both in the case of blunt and sharp noses, and then how to estimate the heat at the wall accordingly to the attitude of the CubeSat.
Generally, as suggested by Koppenwallner et al. [13], when the Mach number is not hypersonic, say < Ma 6, the melting process fin- ishes and the structure starts to cool down and so Q ExchAir is set to zero.
Thermodynamic characterization of the blunt nose flowfield.The hypothesis of local thermodynamic and chemical equilibrium implies that the characteristic time for the chemical reactions and/or vibrational energy to approach equilibrium is significantly smaller than the characteristic flow time.Thus, the flowfield is able to immediately adapt its local chemical composition to the imposed local thermodynamic conditions.The thermodynamic properties of high temperature equilibrium air can be tabulated in form of polynomial correlations as given by Tannehill and Muggie [36] in order to speed up their evaluation and avoid time consuming computational procedures.In particular in this work, we are interested in the following set of correlations: • Temperature as function of density and pressure: = T T P ( , ) Fig. 4. Thermodynamic properties of high temperature air in chemical equilibrium as tabulated by Tannehill and Muggie [36].computed as • specific enthalpy as function of density and pressure: = h h P ( , ) as where ˜is correlated as: • entropy as function of internal energy and density: = s s e ( , ) as 10 can be found in Refs.[16,36,37].Whenever an inverse relation is necessary, the above correlations are numerically solved through a non linear root solver.The three tables relative to the three selected correlations are also plotted in Fig. 4.
Once the high temperature air properties are available, it is possible to compute the thermodynamic conditions across the shock wave.In front of a blunt nose, a bow shock wave detached from the body is present, which can be approximated as normal in the neighbourhood of the stagnation point.Therefore, the iterative procedure suggested by Anderson [16] for the normal shock wave is used.From the balance of mass, momentum and energy across the shock wave, the following set of three non linear equations is obtained: (29) where the subscript means upstream and 2 means downstream the shock wave.All the upstream conditions are assumed to be known because they are provided by the selected atmospheric model whereas the downstream conditions are the unknowns.A fixed-point [38] based iterative procedure is set up, considering the ratio / 2 as the updating variable.The initial guess is computed from the perfect gas relation with specific heat ratio = 1.4: From the momentum balance in eq. ( 30) and the energy balance in eq. ( 31), pressure and enthalpy are computed, respectively.The obtained values are used to invert eq. ( 26) and get thus the density 2 .Dividing it by , the updated value of the ratio / 2 is computed and used as input for the successive iteration.The procedure converges in few iterations accordingly to the selected tolerance.Eq. ( 25) provides the values of the temperature T 2 .The temperature jump across the shock is plotted in Fig. 6 for different values of upstream pressure.Finally, the internal energy e 2 is computed through the First Law of Thermodynamics: and eq.( 28) provides the entropy value s 2 .The thermodynamic conditions at the stagnation point, indicated with the subscript t2, are computed following Bertin [26] and with reference to Fig. 5 for the nomenclature.The total enthalpy can be computed equivalently from the upstream or downstream conditions: and, since the transformation is isoentropic, it follows that = s s t2 2 .Then, the following set of non linear equations can be solved numerically:    (37) to compute the three unknowns: e t2 , P t2 , t2 .Applying again eq. ( 25), it is possible to compute the total temperature at the stagnation point T t2 .The dynamic viscosity is computed through the Sutherland law: where the coefficient has been added to take into account the high temperature of the flow at the edge of the bounary layer [26].Since the static pressure is constant across the thin boundary layer, we can set , , where now the subscript w indicates the conditions at the wall.Eq. ( 25), set as: can be inverted to compute the density at the wall w t , and eq.( 26) to get the enthalpy at the wall: The Sutherland law [26] provides the dynamic viscosity at the wall: In addition, in order to estimate the dissociation enthalpy at the edge of the boundary layer, the composition of the gas is required.A very simplified iterative process is here set up as described by De Luca [39] and by Anderson [16].A mixture of thermally perfect gases is considered including O, O 2 , N, N 2 and NO.Thus, given the working temperature T t2 and pressure P t2 , the unknowns are the number of moles (per unit volume) of these 5 species: n O , n O2 , n N , n N2 , n NO , according to the following chemical reaction: where the symbol N indicates the gram-atom number per unit volume of a specific species initially available.These quantities must be available for the geographical position and altitude of the spacecraft.For instance, the atmospheric model NRLMSISE-00 (sec.5.1) selected in this work provides such type of information.Since only two atomic species are present in the mixture: oxygen and nitrogen, two mass conservation conditions in terms of gram-atoms can be written: Since the unknowns are 5 and so far only 2 equations are given, the 3 remaining equations of chemical type can be provided by the standard formation equilibria: where i indicates the total variation (products minus reactants) of the stoichiometric number of moles associated to the specific i-th reaction: , .The equations are written exploiting the relation between the constant of standard formation in terms of partial pressures K P and the constant of standard formation in terms of number of moles K n : K P and K n can be specified for the three selected reactions as: Under the hypothesis of chemical equilibrium, K P is a function of temperature only.It can be computed from statistical thermodynamics or can be included in terms of tables.The JANAF tables are a suitable source, accessible from Ref. [40].Therefore the only additional unknown that has been introduced is the total number of moles n tot .Its definition closes up the balance between unknowns and equations: To solve the non linear system, composed by equations ( 43) and ( 44) and ( 49)-( 51), the fixed-point based iterative procedure suggested by De Luca [39] has been implemented.The initial guess for n tot is 1.1.Actually, this non linear system can be reduced by substitution to a system of only 2 equations, the mass conservation equations 43 and 44, and 2 unknowns, n O2 and n N2 , which is solved numerically.Finally, the updated value of n tot is provided by its definition.The procedure converges few iterations accordingly to the selected tolerance.The composition at = P atm 1 is shown in Fig. 7. Once the solution is reached, the average properties of the mixture can be computed.It is also possible to add other chemical species that have not taken part to the dissociation process but that are reasonably present in air.Indeed, NRLMSISE-00 also outputs the concentration of Argon, Helium and Monoatomic Hydrogen, for a total of = N ˜8 chemical species and a total number of moles He .So the molar fractions X i Fig. 7. Dry air composition at = P atm 1 computed with the iterative procedure described in this report and data from the JANAF tables [40].are: and the average molecular weight : where i is the molecular weight of the i-th species.Finally the mass fractions c i are computed as: Thermodynamic characterization of the sharp nose flowfield.When the spacecraft is pointing the flowfield with an edge or a corner, the nose cannot be reasonably considered blunt.So it is more likely that a shock wave forms attached to the nose tip, oblique or normal accordingly to the deflection angle imposed by the body and to the Mach number.In case a solution does not exist for the straight oblique shock case, the same iterative process introduced in the previous paragraph for normal shock wave is used to characterize the thermodynamics conditions across the shock.When a solution for the oblique shock exists, it can be computed numerically.Here, we suggest an approach based on two nested root-finder loops.
The internal loop is a fixed-point scheme, used to compute the thermodynamic variables across the shock, assuming to know the shock deflection angle i .It is set identically to the case of the normal shock waves, i.e. considering / 2 as updating variable, but substituting the velocity modules V and V 2 with the respective components normal to the shock wave: . i is the deflection angle imposed by the i-th panel, computed as in eq. ( 7), and i is the relative shock deflection angle.See Fig. 8 for clarification.The polynomial correlations given by Tannehill and Muggie [36] to account for the high temperature air properties are used also in this case.
The external loop provides the internal loop with the value of the shock deflection angle i .Despite the hypothesis of local thermodynamic and chemical equilibrium air, the velocity component tangential to the shock wave preserves across the shock [16]: In light of this, the θ-β-V relation, given in Anderson [16], can be derived: which holds true when an oblique shock solution exists.So, let us consider the following function of ˆ: which is zero when ˆis equal to the actual deflection angle i .When f ( ˆ) has no roots in the range [ , 90 ]   i o , an oblique shock wave solution does not exist and so we directly rely on the normal shock algorithm.When the oblique shock solution exists, f ( ˆ) has two roots, relative to the wake-shock and strong-shock solutions, and the (only) maximum of f ( ˆ) is positive.We are interested in the wake-shock solution [16], which is the one corresponding to the smaller value of i .Increasingly progressively ˆfrom = ˆi to = 90 o , f ( ˆ) is evaluated at each step k.
, then the wake-shock solution exists and lies between ˆk and + ˆk 1 .So they can be used either as limits for the initial interval of a bisection root finder method or as initial guess for a Newton-based root finder method [38].Otherwise, if 1 , then the oblique shock wave solution does not exist because the maximum of f ( ˆ) has been crossed without passing through the axis = f ( ˆ) 0, so we rely on the normal shock wave algo- rithm.Based on some experience with this procedure, a 2 or 3 degrees step length is a good compromise between number of f ( ˆ) evaluations and goodness of ˆk as initial guess of the root finder.Some results computed with this double loop root-finder scheme are provided in Fig. 9, varying both the imposed deflection angle and the upstream velocity.
In both cases of edge and corner pointing, opportunely tuned correlations developed for the simple case of the flat plate can be used to estimate the heat at the wall.However, differently from the stagnation point case, a simple closed form correlation is not available under the hypothesis of equilibrium chemical reacting gas.In order to avoid very time-consuming procedures to characterize the boundary layer, we rely on the assumption of thermally perfect gas and consider the flowfield solution for laminar boundary layer described by Eckert [25] and recalled by Dobarco-Otero et al. [11].Indeed, as indicated by other investigators (see for example Moore [41]), a certain percentage of dissociation does not affect the heat transfer at the wall if the dissociation does not occur in the flowfield layers immediately close to the surface.Whence, the Eckert solution can be fully regarded within the target accuracy of this work [11].The engineering idea that strongly helps in simplifying the problem is to consider a reference temperature condition to evaluate the properties of the fluid.This allows to use straightforward analytic relations, that in some sense remind the selfsimilar solutions found for the simpler case of calorically perfect gas (see for instance Anderson [16]).Since the specific heat of the flow may vary within the boundary layer, we compute the reference condition, indicated with the superscript *, through the reference enthalpy as: where the subscript e is used to indicate the conditions at the edge of the boundary layer, approximated as the conditions downstream the shock wave, that is = h h e 2 and = V V e 2 .In addition, since the boundary layer is very thin in the hypersonic regime, the pressure is considered constant and equal to downstream the shock wave: So the density at the wall is computed inverting eq. ( 25), set as: and eq.( 26) gives the enthalpy at the wall: where the Prandtl number is introduced at the reference temperature T * .Experimental data of Prandtl number varying with temperature are taken from Van Driest [17] and reported in Fig. 10.
The pressure downstream the shock wave is also taken as reference pressure: P P 2 * .Thus, eq. ( 26), set as: is inverted to compute the reference density ρ*, while eq.( 25) gives the  to account for the dissociation effects.Heat transferred to the structure as function of the attitude mode.In order to estimate the heat at the wall, which is significantly affected by the bluntness/sharpness of the spacecraft nose, a suitable correlation must be selected accordingly to the specific attitude of the CubeSat.The face pointing mode has a flat nose, thus very blunt.The Fay and Riddel [24] correlation for equilibrium boundary layer is appropriate in this case if slightly modified.It expresses the convective heat power transferred to the structure at the spherical stagnation point under the hypothesis of laminar boundary layer as: where Pr w is the Prandtl number at the wall temperature T w and Le is the Lewis number taken equal to 1.4.The average dissociation enthalpy at the edge of the boundary layer h D is computed from the composition of the gas: where c i are the mass fraction of the species at the edge of the boundary layer and h ( ) o is the standard enthalpy of formation of the species per unit mass, available in [40].The term u x (d /d ) e t2 is the velocity gradient at the stagnation point.Usually a suitable approximation is obtained considering inviscid flow and applying the Newtonian theory to get: The parameter R eff is the effective radius of curvature of the nose.An effective radius has been introduced because an actual curvature radius does not exist for a cubical shape in face pointing mode.A suitable correction coefficient is provided by the Boison and Curtiss's correlation [42] in Fig. 11  .Considering the worst case of maximum heat into the structure, the minimum between the two half dimensions of the CubeSat face is selected as r * .Therefore, at the end: modifying H and W according to the attitude configuration.The contribution due to the radiation with the air flow is generally not very significant for a re-entry from LEO [16].However, it is here considered and computed using the correlation suggested by Hamilton et al. [43]: ) 0.3542 0.5646 (0.306 0.066 )log ( ) where V is introduced in km s / , in kg m / 3 , R eff in meters and q s RadAir is calculated in W cm / 2 .Once the stagnation heat is computed, the following formula, introduced by Koppenwallner et al. [13], can be used to estimate the heat power over the wetted ( > 0 i ) surfaces of the CubeSat: The edge pointing mode has a sharp nose, so the Eckert solution [25] is suitable for this case, assuming a constant pressure along the body surface, two dimensional flow (wedge) and isothermal wall.Thus, we can express the heat power at a given distance x from the leading edge of the i-th panel as: where the coefficient a i is given by: The heat power profile can be integrated over the specific surface of the CubeSat to get the average heating rate q w Conv i over the i-th panel.The averaging formula depends on the orientation of the panel.Let us consider for example that the W edge points to the flow.If the panel is adjacent to W, the average heat power is given by: If the panel is perpendicular to the pointing edge, the inclination angle i is zero.The average heat power is calculated as: where the integral is approximated numerically [38].Note in particular that in the two equations above the three dimensions W, H and L of the CubeSat have been introduced for the case having the W edge pointing to the flow.When the other two edges are taken pointing to the flow, the formulas have to be modified accordingly.The heat power due to the radiation with the air flow is added using again the Hamilton's correlation and the Koppenwallner's formula: = + + q q q (0.1 0.9 sin ) The corner pointing mode is computed similarly to the edge pointing mode.The difference is that now the flow cannot be approximated as two dimensional and a coefficient to account for the three-dimensional relieving effect of the flow must be added.Dobarco-Otero et al. [11], but also Anderson [16] and Bertin suggest to multiply by 3 the a coefficient for the two dimensional flow above the wedge, to obtain the a coefficient for the three dimensional flow above the cone.The same approach is considered here, as well.Thus, for the corner pointing we use: The average heat power is integrated on the three wetted surfaces using eq.( 74), where the dimensions are selected according to the specific panel under consideration.Eq. ( 75) is then used to take into account the radiation contribution.Finally, for all the pointing modes described so far, a discrete summation on the facets of the CubeSat, as done in the Free molecular regime in (23), gives the total heat power Q ExchAir entering into the body due to the thermal exchange with the air flow.
The tumbling mode is approximated with the same averaging approach used in the aerodynamic model for the drag coefficient.More specifically, the total heat power for the tumbling mode is computed averaging among the other and weighting on the number of faces, edges and corners: (77)

Dynamic model
Following the usual approach of an object oriented code [15,44], the spacecraft is considered as a point mass, i.e. three degrees of freedom, flying inside the rotating terrestrial atmosphere.The attitude mode depends on the specific spacecraft.If particular components or appendages of the spacecraft determine a certain stabilization (as when the spacecraft is equipped with the D3 system, see section 7 for details), a specific attitude mode can be selected.Otherwise, the spacecraft can be considered randomly tumbling, modeled as weighted average among all the other modes, as explained in the previous sections.Thus, the information on the local angular velocity is not required and the integration of the three moment equations (Euler equations [19]) can be avoided.Therefore, the numerical simulation is described by only eight first order differential equations: three kinematic equations, three dynamics equations, an equation determining the mass variation in time due to the melting process and a thermal equation to account for the temperature increasing of the body.These equations are presented in the next paragraphs following the work of Vinh at al [18].

Kinematic and dynamics equations
Firstly, with reference to Fig. 12, let us define: • the local horizontal plane as the plane passing through the vehicle and perpendicular to the position vector r • the local vertical plane as the plane containing the position vector r and the velocity vector v relative to the Earth's rotating atmosphere.
The equations of motion can be written in a suitable scalar form for the re-entry dynamics if they are written in the Vehicle-Centred Intrinsic (VCI) reference frame.Defined by the three Cartesian unit vectors i j k { , , }, it has: the origin in the vehicle center of mass, the second axis oriented as the relative velocity vector v, the first axis in the local vertical plane perpendicular to v and the third axis in the local horizontal plane completing the right-hand rule.In this frame, we can write [18]: where: • r and v are the magnitude of the position and of the relative velocity vectors, respectively; λ and φ are longitude and geocentric latitude; • γ is the flight path angle defined as the angle in the local vertical plane between the local horizontal plane and the relative velocity vector v, positively when v is above the local horizontal plane; • ψ is the heading angle defined as the angle in the local horizontal plane between the local parallel (i.e. the unit vector j ) and the projection of the relative velocity vector v in the local horizontal plane, positively in the right-handed direction around the position vector r.
• m is the overall mass of the spacecraft; • = × rad s 7.292115 10 / 5 is the angular velocity of the Earth rotation.
H are the three components of the overall force vector acting on the vehicle, expressed in the VCI reference frame.force vector F is the sum of two contributions: • the aerodynamic drag: where is the atmospheric density, which is described in this work by the 2001 United States Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar Exosphere model (NRLMSISE-00) [45]; is the free stream velocity that coincides with the magnitude of the relative velocity vector: V v since the presence of wind is totally β is the ballistic coefficient, defined as: where c D is the drag coefficient given by the aerodynamic model (equations (( 6), ( 8), ( 15) and ( 18) and Fig. 3 as function of Knudsen number, Mach number and attitude mode) and = S HW Ref is the reference area used to define the drag coefficient, in this case the frontal area of the CubeSat.
• the gravitational force: (86) which includes the effect of the second zonal harmonic J 2 , expressed as given by Vallado [46].Specifically, µ is the Earth gravitational parameter, = J 0.0010826269 2 is the second zonal harmonic coefficient, i is the orbital inclination and u is the argument of latitude (sum of argument of perigee and true anomaly).The vectors triad R { ˆ, , ˆ} identifies the Earth-Centred Orbital reference frame (ECO).The initial condition for all the simulations in this paper is defined De-Orbit (DEO) point.It is placed on a 120 km circular orbit, having inclination of 45°, Right Ascension of the Ascending Node (RAAN) of 0°, argument of latitude of 180°.Thus, the initial latitude is above the equator and the initial longitude is set to 99°to have a re-entry over the South Pacific Ocean.

Thermal and melting equations
As a worst case, we consider that the internal payload does not play a role in the thermal balance, i.e. it is approximated as thermally isolated, and for simplicity the external Aluminum case of the CubeSat is considered as a single thermal node.Thus we can write that: where the total heat power Q is provided by the aerothermodynamic model (equations ( 21) and ( 77)).The other variables:

and c w
Al are the melting temperature and the specific heat of the Alluminium wall.c w Al is considered temperature dependent as given in Lips and Fritsche [15] for ORSAT: given in J kgK /( ).m th is the thermal mass, that is only the mass of the external case, calculated as: where δ is the thickness of the structure, which decreases with time because of the melting process, = kg m 2700 / w Al 3 is the density of the Aluminum wall and W, H and L are the CubeSat dimensions, considered constant.The thickness δ is an independent variable of the problem computed through the integration of the melting equation, as given by Koppenwallner et al. [13]: where is the heat of ablation of Aluminum and S Ext is the external surface of the spacecraft, calculated as in eq.(20).Numerically integrating in time the 8 ordinary differential equations: (78) to (83) together with (102)-( 103) and ( 105)-(106), it is possible to fully characterize in an approximate way the re-entry of a CubeSat flying with a particular pointing attitude or randomly tumbling.The integration is stopped when the condition = 0 is reached which implies that the thermal mass of the spacecraft is completely melted.

Results for CubeSat re-entry
In Fig. 13 the variation of the most significant quantities: geodetic altitude, relative velocity modulus, drag coefficient, net heat power rate, temperature profile, is plotted along the integration of the re-entry dynamics.In particular, a 1U Cubesat ( = W cm 10 , = H cm 10 , = L cm 11.3 ) is considered and all the attitude modes are analyzed.In order to show all the velocity regimes up to ground impact, the melting process is numerically forced to be negligible setting the heat of ablation very high (h Al

Abl
).For what concerns the drag coefficient (see Fig. 13c), we can recognize an opposite trend between Free molecular and Continuum regimes.In the former, the shear stress provided by the Schaaf and Chambre's model represents good percentage of the total drag force and it is particularly relevant when the spacecraft is in edge or corner pointing.In the latter, the shear contribution is neglected by the Newton theory, so the highest drag is provided for the face pointing mode.Since the atmospheric impact is particularly sensitive to the drag experienced at high altitudes, the CubeSat will be slowed down by the atmospheric drag a bit sooner if it flies in corner or pointing with respect to face pointing mode (Fig. 13b).Consequently, it reaches sooner the ground (Fig. 13a).
The net heat power (convection minus re-radiation) is minimum in face pointing mode because it is always inversely proportional to the bluntness of the nose (Fig. 13d).Corner and edge pointing are sharp noses and so they get much more heat.Whence the spacecraft reaches later the melting temperature flying in face pointing mode than flying in the other modes (Fig. 13e).The modeling of the tumbling mode gives effectively a good average among the other modes.In this case, the tumbling mode closely approaches the edge pointing mode.
In Fig. 14 a larger CubeSat (6U) is considered, having all the dimensions unequal (3.2.1 configuration in Table 2 and = L cm 34 .Thus, it is possible to notice how the drag coefficient (Fig. 14a, c and e) and the net heat power (Fig. 14b, d and f) change accordingly to which specific face or edge is selected.For this type of CubeSat, the largest drag coefficient is experienced flying with the WL face pointing to the flow (Fig. 14a).This is an expected result because the WL face pointing mode has the largest frontal area.The heat, instead, is minimum in this mode because the effective radius of curvature is the largest (Fig. 14b).
Another comment can be made regarding the small discontinuities visible in Fig. 14d as well as in Fig. 14f.They occur because, as long as the Mach number decreases, an oblique shock wave solution may not exist any more.So the code switches to the normal shock wave algorithm considering a detached bow shock wave.This is in some sense equivalent to a larger bluntness of the nose and so it implies a reduction of the heat getting into the body.
In Fig. 15 the melting process is, instead, considered active and all the CubeSat configurations in Table 2 are analyzed.When the temperature of the spacecraft reaches the melting one, the thermal mass starts decreasing with time.When it becomes zero, the integration stops and the final geodetic altitude reached by the spacecraft is the demise altitude h Dem .
In Fig. 15, we see that the demise altitude is always around km 93 98 and we can recognize two trends: 1. Increasing the length of the CubeSat and consequently the mass, the average ballistic coefficient increases.So the re-entry is shallower and the CubeSat demise occurs later in time but roughly at the same altitude; 2. Increasing the surface to frontal area ratio, the CubeSat experiences on average more heat especially in the Free molecular regime.Consequently, it demises at higher altitudes.
In Fig. 16, we compare few test cases with some results of ORSAT and SCARAB simulations taken from Ref. [29] of boxes re-entering from a km 122 altitude on a slightly elliptical orbit.Note that in Ref. [29] several test cases of spheres and cylinders are analyzed and compared.Since these shapes do not directly pertain with this work, we focused on the boxes shape only.Unfortunately, the plots of only four different boxes are published in Ref. [29], so we can base the comparison only on the available data.The precise initial conditions of the simulations as well as the geometrical dimensions and the physical properties of the boxes are reported in the referenced paper.A significant discrepancy in time was already recognized between ORSAT and SCARAB results [29], clearly visible in Fig. 16a and in Fig. 16c.The same time effect is also evident in the results computed with the model presented in this work.Nevertheless, when the altitude is plotted versus the velocity modulus (see Fig. 16b) it is found a good agreement.This implies that the deceleration rate is almost in accordance.In addition, the simplified thermal model with a single thermal node adopted in this work does not allow to estimate the maximum temperature on the box structure, accounting only for an average value.However, the average wall temperature estimated by the model closely approaches the maximum temperature computed by ORSAT and SCARAB when the temperature reaches almost the melting one (see Fig. 16c and d).This is expected since at that point the temperature on the body is very high and almost averaged close to the melting one.This implies a good agreement in the estimated thermal loads, as well.

Model of the D3 system
When the CubeSat is equipped with the D3 system, its attitude is 3axis stabilized with the WH face pointing the flow [9].Thus, the initial part of the re-entry trajectory is always accounted having HW face pointing.The bending stiffness of the D3 booms has been designed in  with fully deployed (about m 3.75 ) booms.The torque experienced by the D3 booms is computed along the numerical simulation, checking when it exceeds C Max .When this condition is verified, the D3 booms are considered fully bent and parallel to the flow, contributing to the drag force only through the shear stress with the surrounding flow.The numerical simulation proceeds till the first melting condition is reached.Indeed, two possible situations may occur: • if the spacecraft structure melts first, the integration stops and the final reached altitude is the demise altitude h Dem ; • if the D3 booms melt first, the integration proceeds considering the CubeSat only, in tumbling mode; then, when also the spacecraft structure melts the integration stops and the final altitude is the demise altitude.
All these conditions are illustrated in Fig. 17.In the following paragraphs, we briefly describe how to model the D3 booms and how to compute the experienced torque and the heat power that enters into the D3 booms mass.All the selected correlations are collected in Table 3.

Aerodynamic model
The four D3 Booms are modeled as four additional panels of the spacecraft, having only one face pointing to the flow and inclined with respect to it of a 70deg angle.all the correlations used for the CubeSat are implemented for the D3 booms as well.In particular, in the Free molecular regime the Schaaf and Chambre's analytic model [21] is used with = = 0.9 . Thus, the drag coefficient for the single D3 boom is calculated as: where the superscript D3B refers to the D3 boom.In the Transition regime the Wilmoth bridging formula [30] is used.In the Continuum regime, a distinction must be made on whether the D3 booms are bent or not.When they are still deployed, the Modified Newton Law [22] gives a good approximation of the pressure on the booms and the shear stress is neglected.On the contrary, once they are bent, the drag is a result entirely of the shear.To account for it, the strong interaction theory of Li and Nagamatsu introduced in Ref. [47] and recalled by Boettcher et al. [48] is suggested.Integrating the local skin friction coefficient, the drag coefficient for the single boom can be expressed as given by Koppenwallner et al. [49]: where c T T ( / ) f w t 0 2 is a function of the wall-to-stagnation temperature ratio and of the specific heat ratio.It is tabulated in Ref. [47] for = 1.4.C is the Chapman-Rubesin constant computed as: Re is the free stream Reynolds number based on the D3 boom length: Finally the overall drag coefficient is computed as: where the superscript Cub refer to the CubeSat.The reference area S Ref is always the frontal surface of the CubeSat.The torque C D3B experienced during the re-entry by the D3 boom is finally computed as: (98)

Aerothermodynamic model
The four D3 booms are modeled as a single thermal node, isolated from structural body.The choice of using a single thermal node is due to the fact that the entire boom is subjected to the same airflow.Therefore, the heat input depends only on inclination angle with respect to the airflow.Since the booms are approximated as perfectly straight from root to tip, the heat input is constant along the boom length and the temperature increases mostly 3 uniformly.The hypothesis of isolation from the CubeSat is introduced because the thickness of the booms is very small (about mm 0.0762 ) and consequently the heat transferred by conduction to the spacecraft is negligible.The heat power is computed using the same correlations used for the CubeSat structural mass.In particular, in the Free molecular regime the Schaaf and Chambre's analytical model [21] is used with = a 1 c .In the Transition regime, the Wilmoth's bridging formula [30] is preferred to the Legge's formula [23] because it foresees a smoother net heat power on the D3 booms.So we write: with ϕ as function of the Kudsen number computed as in eq. ( 9).In the Continuum hypersonic regime, the Koppenwallner's formula [13] (eq. ( 70)) allows to compute the heat power at the wall from the stagnation heat, given by the Fay and Riddel's correlation [24] (eq.( 65)) since, as far as the D3 system is not melted, the CubeSat flies in HW face pointing mode.The Hamilton correlation [43] (eq.( 69)) gives the contribution due to radiation with the air flow.At low Mach number ( < Ma 6), the heat entering into the structure is set to zero.Finally, the total thermal power that enters into the single D3 boom due to the thermal exchange with the air flow is: and the net heat is: accounting for the heat re-radiated towards the environment.= ss is the stainless steel emissivity and T w D3B is the D3 booms wall temperature.

Dynamic model
Adding the D3 booms to the model, two ordinary differential equations must be added to the model, for a total of 10 equations, to account for the temperature increase of the D3 booms and the resulting melting process.Following the same approach used for the CubeSat, we can write the thermal equations as:

Results for CubeSats equipped with the D3 system
The main effect of having the CubeSat equipped with the D3 system is to increase the drag coefficient and consequently to decrease the ballistic coefficient.To verify this, compare for instance the drag profile in Fig. 18c with Fig. 13c.Thus, the longer the booms, the shorter is the time spent to cross the altitude range 120-110 km (Fig. 18a and b).This implies that the deceleration is faster with higher localized peaks of heat, especially at the stagnation point.Nevertheless, the total heat input is lower because the CubeSat spends less time inside the atmosphere and because it flies in stabilized face pointing mode.Consequently, we can expect that the D3 system shall let the spacecraft reach lower demise altitudes, even lower than km 90 with D3 booms deployed for at least m 0.5 , as shown in Fig. 20.Actually, since the booms bending altitude increases with the increase of the booms length, the jump in the drag coefficient occurs at higher altitudes.This is a balancing effect that makes the final demise altitude of the CubeSat have a very slight dependency on the initial length of the booms.In addition, analysing more Fig.20, we can see that if the booms length is smaller than m 0.3 , the melting temperature is reached sooner than the torque force becomes high enough to bend the booms.Thus, the booms melt even before than the CubeSat structure.If it happened, the CubeSat would start tumbling and consequently getting more heat, as shown in Fig. 19a for the case with = L m 0.3 D3B .Therefore, the CubeSat demise altitude slightly increases (Fig. 20).On the contrary, for longer booms length, at least m 0.4 , the heat power per unit of area that gets into the CubeSat structure is much higher than that into the D3 booms.Compare Fig. 19b with Fig. 19a.This is mainly due to the reradiating cooling power of the booms, which occurs upon a double area (front and back) with respect to the convective heating (only front).In addition, the booms are made of stainless steel which can reach very high temperatures before melting.Since the re-radiating cooling increases with the fourth power of the temperature, when the booms temperature approaches the melting one, the re-radiation is able to dissipate a large amount of heat and quickly re-balance the convective heating (Fig. 19b).Whence, it seems very unlikely that, if the booms are long enough, they are going to melt before the CubeSat structure.The model is not able to foresee what happens to the D3 booms, whether they demise or not, after the CubeSat melting event.In Fig. 19d, we can clearly see that as the length of the D3 booms increases, the bending altitude increases as well, because the momentum at the root of the boom imposed by the drag force is larger.When the booms bend, the heat that they experience suddenly decreases because they are considered aligned with the flow and the reradiation with the environment cool them down.Since the thermal mass of the booms is very small, the balance between convective heating and re-radiating cooling is then reached again (Fig. 19b) and the booms temperature reaches very quickly the steady state temperature.This is why, after the bending event, the booms temperature is roughly the same no matter which is their length.The thermal mass of the CubeSat is, Fig. 18.Dynamic-related quantities of the simulation of a 1U CubeSat from a 120 km altitude circular orbit equipped with the D3 system.Five different booms lengths are evaluated with active melting process.instead, higher and it might not have enough time to adapt to the steady state temperature, especially at very high altitudes.Nevertheless, the balancing effect given by the booms bending damps this effect and provides an altitude of demise of the CubeSat only slightly dependent on the boom length (Fig. 20).
Finally, we can evaluate the demise altitudes for different types of CubeSat but keeping constant the initial length of the D3 booms, as illustrated in Fig. 21.The same trends that we recognized for the case without D3 system repeat.The longer the CubeSat, later and at lower altitudes it melts, even though the demise altitude is affected very slightly.A relatively large increase of the demise altitude is, instead, observed with the increase of the surface to frontal area ratio and the CubeSat mass.Indeed, approaching more and more a cubic shape, the heat going into the body is on the average higher.The larger mass implies a larger ballistic coefficient that lets the spacecraft take more time inside the atmosphere and so collect more heat.

Conclusions
The aerodynamic and aerothermodynamic analysis of a spacecraft re-entering into the atmosphere presents several challenges from the modeling standpoint.The simulation has to be divided in ranges as function of the Knudsen and the Mach numbers.For each range, a specific formulation has to be selected according to the spacecraft Fig. 19.Thermodynamic-related quantities of the simulation of a 1U CubeSat from a 120 km altitude circular orbit equipped with the D3 system.Five different booms lengths are evaluated with active melting process.geometry to estimate both the aerodynamic forces acting on the spacecraft and the heating power entering into the structure.This paper presents an organized and comprehensive set of formulations to analyze in an approximate way the re-entry of smallsats having a rectangular prismatic shape.All the selected correlations, laws and data provide the description of the phenomena with the suitable accuracy that pertains the atmospheric re-entry problem, whereas keeping algebraic form.Whence, it is not required to run intensive computational fluid dynamics based analysis, preserving thus the time for a single simulation of the re-entry in a standard desktop computer within few tens of seconds.The test cases simulations provided in this work prove that the overall model is consistent with the underlying physics of the problem.The reduction of the drag coefficient, evidenced also in other studies, passing from Rarefied to Continuum regime is correctly model.In addition, the opposite trend between the two regimes due to the different influence of the shear contribution in the overall budget of the drag force is evidenced.When the spacecraft flies in a sharp nose configuration, the model predicts higher heat input than in more blunt nose configurations, in accordance with the theory.The model is also able to capture the effect of the ballistic coefficient related to different configurations and mass of the CubeSat.Increasing the CubeSat length, the demise occurs later but the altitude in not substantially affected.On the contrary, the increase of the surface to frontal area ratio, increases the demise altitude.Adding the D3 system has the effect to strongly increase the drag experienced by the spacecraft and so to shorten the reentry and make it steeper.Furthermore, the D3 system aerodynamically stabilize the spacecraft in face pointing mode.This reduces the heat input into the spacecraft and so decreases the demises altitude lower than km 90 with the D3 booms deployed for at least m 0.4 .The model is also able to approximately estimate at which altitude the booms bend and/or melt, when it occurs before the spacecraft structure melting event.It is shown that whether the booms first bend or melt depends on their length.For the 1U CubeSat, it is probable that they melt before bending if the booms length is lower than m 0.3 .Nevertheless, the booms length does not significantly affect the CubeSat demise altitude because of the balancing effect of the booms bending.Finally, we can recognize the same effects due to the change of the ballistic coefficient and to the surface to frontal area ratio that have been highlighted simulating only the CubeSat.The first limitation of the model is that it is not able to estimate the heat flux that gets into the CubeSat payload, because of the simplified thermal model.This estimation is outside the scope of this work but it can be regarded as an interesting further development.The second limitation concerns the demise altitude of the D3 booms when it occurs after the CubeSat demise event.Since the attitude of the D3 boom flying alone is not known a priori, the approximated models described in this work do not allow one to estimate its aero-and aerothermo-dynamic parameters.Future work includes the update of the models to take into account this scenario.
) where c D FM and c D CR are the drag coefficients computed considering the flowfield in Free Molecular regime or in Continuum regime, respectively.The parameter ϕ is given as function of the local Knudsen number by: coefficients a 1 and a 2 should be adjusted to give the best overall description of the transitional flow.In order to get a continuous and smooth behaviour, they are taken: molecular-Transition and the Transition-Continuum Knudsen boundaries.
c D HW , c D WL , c D HL accordingly to which specific face points the flow; for edge pointing among: c D W , c D H , c D L accordingly to which specific edge points the flow and for corner pointing only one coefficient is necessary: c D C

Fig. 3 .
Fig. 3. Drag coefficient as function of upstream Mach number for rotating cubes fired through a ballistic range.Taken from Fig. 17 chapter 16 of [35].

Fig. 5 .
Fig. 5. Illustration of the shock wave in front of the CubeSat in face pointing mode (blunt nose).The subscripts nomenclature used in the thermodynamic analysis of the flow is highlighted.

Fig. 6 .
Fig. 6.Temperature jump across a normal shock wave versus upstream flow velocity and pressure with = T K 225 and considering high temperature air in chemical equilibrium.

Fig. 8 .
Fig. 8. Illustration of the oblique shock waves in front of the CubeSat in edge pointing mode (sharp nose).The subscripts nomenclature used in the thermodynamic analysis of the flow as well as the oblique shock wave geometry are highlighted.

Fig. 9 .
Fig. 9. Properties across an oblique shock wave as function of the upstream flow velocity, with = P atm 0.01 and = T K 225 , and considering high temperature air in chemical equilibrium.
that accounts for the truncated flat nose.Since for the flat face =

Fig. 12 .
Fig. 12. Reference frames illustrations: Earth Centred Earth Fixed (ECEF) frame in red, Earth Centred Local (ECL) frame in blue and Vehicle centred Intrinsic (VCI) frame in green.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 13 .
Fig. 13.Simulation of a 1U CubeSat from a 120 km altitude circular orbit considering infinite heat of ablation (h Al

Fig. 14 .
Fig. 14.Simulation of a 6U 3.2.1 CubeSat from a 120 km altitude circular orbit considering infinite heat of ablation (h Al Abl

Fig. 15 .
Fig. 15.Simulation of eight different types of CubeSats in tumbling mode and active melting process.

Fig. 16 .
Fig. 16.Comparison with ORSAT and SCARAB of few test cases simulation of boxes re-entering from a km 122 slightly ellitptical orbit.ORSAT and SCARAB data are taken from Ref. [29].

3
and the D3 booms thickness D3B is computed integrating the melting equation:

Fig. 21 .===== 2 =
Fig. 21.Simulation of eight different types of equipped with the D3 system having m 0.5 booms length and active melting process.

Table 1
Summary of all the correlations used in the aerodynamic and aerothermodynamic models.

Table 3
Summary of all the correlations used in the aerodynamic and aerothermodynamic models of the D3 system. ]