Atmospheric Reentry Dynamics of Conic Objects

Abstract

One of the key issues in a reentry risk analysis is the calculation of the aerodynamic coefficients. This paper presents a methodology to obtain these coefficients and couple it to a code that computes re-entry trajectories considering six degrees of freedom. To evaluate the different flight conditions encountered during the natural re-entry of conical objects, the Euler Equations for gasdynamics flows are used. A new scheme TVD (Total Variation Diminishing) is incorporated to a finite volume unstructured cell-centred formulation, for application to three-dimensional Euler flows. Finally, numerical results are obtained for a conical body at different attack angles and Mach. With these results, the calculation of the trajectories during atmospheric re-entry is completed.

1. Introduction

In the case of natural reentries (non-controlled), the orbital evolution of an object can only be monitored, with no or limited ability to control risks. The time window for reentry of a satellite is usually provided with a standard error of ±10% to ±20% of the remaining orbital lifetime. For the controlled reentries, it is required to simulate the different scenarios until the right window for the mission is found, being that the total or partial disintegration, or the landing on a safe place.

In this paper, the main objective is to conduct numerical simulations of the supersonic flow regime on a conical body, thereby using a code developed at the Department of Aeronautics of the National University of Córdoba, Argentina [1, 2]. This code uses the technique of finite volumes for solving Euler equations. The spatial discretization of the domain is done through a mesh of unstructured tetrahedral volumes. It has implemented a new technique for choosing the limiter functions that can reduce the artificial viscosity without the loss of strength robustness of the Total Variation Diminishing scheme—TVD [3–6].

The drag and lift coefficients, 𝐶_𝐷 and 𝐶_𝐿, obtained by the numerical simulation of compressible flows are used in a code that allows to evaluate the trajectories considering six degrees of freedom [7]. As a result of this research, the trajectories of reentry into the Earth’s atmosphere for conical objects having different initial flight conditions are presented.

2. Methodology

2.1. Description of the Numerical Scheme for Compressible Flow

The three-dimensional Euler equations can be written as

𝜕𝐔𝜕𝑡+∇⋅𝐅=0,(2.1) where U is the vector of conservative variables, and F is the 3D vectorial flow.

The temporal change of the conservative variables can be expressed as

𝐔^𝑛+1=𝐔^𝑛−Δ𝑡𝑉𝑜𝑙^{𝑙faces𝑙=1}𝐅^∗𝑙⋅𝐧_𝑙𝐴_𝑙,(2.2) where the flux of the conservative variables F has been replaced by the numerical flux tensor 𝐅^∗. Vol indicates the volume where the integration is performed, 𝐧_𝑙 is the outward normal to the control surface (𝐴_𝑙).

Equation (2.2) allows the use of a locally aligned system of coordinates whose unit vector i coincides with the normal to the face 𝑙 of the cell, and the unit vectors j and k are tangential directions. To achieve second-order accuracy, the numerical flux at the interface between cells l and l+1 in the direction normal to the face 𝑙 is calculated by [8]

𝐟^∗𝑖+1/2=𝐟^𝑖+𝐟_𝑖+12+12^𝑚Φ^𝑚𝑖+1/2𝐊^𝑚𝑖+1/2,(2.3) where f _i and f _i+1 are the physical fluxes normal to the face in each cells, 𝐊^𝑚𝑖+1/2 is the m-th right eigenvector, and Φ^𝑚𝑖+1/2 is, in the original Harten-Yee scheme [9–11], defined as

Φ^𝑚𝑖+1/2=𝑔^𝑚𝑖+𝑔^𝑚𝑖+1−||𝜆^𝑚𝑖+1/2+𝛽^𝑚𝑖+1/2||𝛼^𝑚𝑖+1/2,(2.4) being

𝑔^𝑚𝑖=𝑆2max0,min||𝜆^𝑚𝑖+1/2𝑆=sign𝜆𝛼^𝑚𝑖+1/2_𝑚𝑖+1/2||,𝑆||𝜆^{𝑚𝑖−1/2}||𝛼^{𝑚𝑖−1/2},(2.5)𝛽,(2.6)_𝑚𝑖+1/2=⎧⎪⎨⎪⎩𝑔^𝑚𝑖+1−𝑔^𝑚𝑖𝛼^𝑚𝑖+1/20if𝛼^𝑚𝑖+1/2if𝛼_𝑚𝑖+1/2≠0=0,(2.7) where 𝛼^𝑚𝑖+1/2 is the jump of the conserved variables across the interfaces between cells i and i+1, 𝜆^{𝑚𝑖−1/2} is the m-th eigenvalue of the Jacobian matrix, 𝑔^𝑚𝑖 is the limiter function, and S is the sign function of the corresponding eigenvalue [8–11]. Since the local Riemann problem is solved with rotated data, the eigensystem is calculated in the locally aligned coordinate frame.

The limiter function given in (2.5) is known as minmod [3–6]. The minmod selects the minimum possible value, so that the scheme is TVD. The other end is the limiter function superbee that ponders the contribution of the high-order flux [3]. The only implementation of the superbee function leads to an excessively compressive scheme which it is not very robust for general practical aerospace applications [6].

In the numerical solution of the three-dimensional Euler equations, five wave families appear. If the five wave families are enumerated in correspondence with their speed, being one the slowest and five the fastest, it can be demonstrated that for waves of the families two to four, the characteristic velocities at both sides of the discontinuity are the same and equal to the velocity discontinuity [3, 5]. This property makes it very difficult to solve theses waves accurately. Generally they are solved diffusely because the numerical methods incorporate a large amount of artificial viscosity to track the contact discontinuity.

In this work, the possibility of implementing different limiter functions for different wave families is explored. The objective is to improve the numerical resolution of the discontinuities associated with the families two to four using compressive limiter functions (superbee), and without losing robustness mainly due to the use of diffusive limiter functions (minmod) for the wave families one and five. This technique implements the utilization of the superbee limiter function only in linear degenerate waves and the minmod function in nondegenerated nonlinear waves [1, 2].

To introduce in the numerical fluxes calculations the limiter function superbee, (2.5) is replaced by the following expression:

𝑔^𝑚𝑖=⎧⎪⎨⎪⎩0if𝛼^𝑚𝑖+1/2max[0,min(2𝑘,1),min(𝑘,2)]1𝛼^{𝑚𝑖−1/2}<02||𝜆^{𝑚𝑖−1/2}||𝛼^{𝑚𝑖−1/2}if𝛼^𝑚𝑖+1/2𝛼^{𝑚𝑖−1/2}≥0(2.8) being

𝑘=||𝜆^𝑚𝑖+1/2||𝛼^𝑚𝑖+1/2||𝜆^{𝑚𝑖−1/2}||𝛼^{𝑚𝑖−1/2}.(2.9)

To improve the overall scheme robustness, the implementation of different limiter functions is carried out only in those cells interfaces where the greater relative intensities of the discontinuities in central waves are registered, and using the conventional Harten-Yee TVD scheme in all other cases. Notice that the comparison among the intensity of the waves cannot be made using directly the coefficients of the spectral decomposition (𝛼^𝑚𝑖+1/2) since these coefficients depend on the module assigned to each eigenvector.

In the local coordinate system adopted for computing the numerical fluxes across each face, the corresponding eigenvectors are given by [5]

𝐊₁=⎡⎢⎢⎢⎢⎢⎣𝐻𝑢𝑤−−𝑣1𝑢𝑐𝑐⎤⎥⎥⎥⎥⎥⎦,𝐊₂=⎡⎢⎢⎢⎢⎢⎢⎣𝑢₂𝑤𝑣1𝑢+𝑣²+𝑤²2⎤⎥⎥⎥⎥⎥⎥⎦,𝐊³=⎡⎢⎢⎢⎢⎢⎣𝑣0010⎤⎥⎥⎥⎥⎥⎦,𝐊₄=⎡⎢⎢⎢⎢⎢⎣𝑤0001⎤⎥⎥⎥⎥⎥⎦,𝐊₅=⎡⎢⎢⎢⎢⎢⎣𝐻𝑢𝑤++𝑣1𝑢𝑐𝑐⎤⎥⎥⎥⎥⎥⎦,(2.10) where 𝐻 is the stagnation enthalpy; u, v and w are the velocity vector components, and c is the sound velocity. It can be deduced from (2.10) that 𝛼^1𝑖+1/2, 𝛼^2𝑖+1/2 and 𝛼^5𝑖+1/2 measure the density jump in the waves 1, 2, and 5, respectively, and that 𝛼^3𝑖+1/2 and 𝛼^4𝑖+1/2 measure the momentum jump in waves three and four. To compare these jumps it became necessary to select reference values for the density and velocity. Thus,

𝐼₁=||𝛼^1𝑖+1/2||𝜌_𝑟𝑒𝑓,𝐼₂=||𝛼^2𝑖+1/2||𝜌_𝑟𝑒𝑓,𝐼₃=||𝛼^3𝑖+1/2||𝜌_𝑟𝑒𝑓𝑢_𝑟𝑒𝑓,𝐼₄=||𝛼^4𝑖+1/2||𝜌_𝑟𝑒𝑓𝑢_𝑟𝑒𝑓,𝐼₅=||𝛼^5𝑖+1/2||𝜌_𝑟𝑒𝑓.(2.11)

In this investigation, 𝜌_ref = 0.5 (𝜌_𝑖 + 𝜌_𝑖+1) is taken as density reference, and as the velocity reference of the average of the sound velocities at the cells 𝑢_𝑟𝑒𝑓 = 0.5 (𝑐_𝑖 + 𝑐_𝑖+1), where 𝑐_𝑖 is the sound velocity. The parameters I_i permit to measure the wave intensities.

Finally, if the maximum of 𝐼₁, 𝐼₅ is higher than the maximum of 𝐼₂, 𝐼₃, 𝐼₄, the conventional Harten-Yee TVD scheme is used; otherwise, the values of 𝑔^2𝑖, 𝑔^3𝑖, 𝑔^4𝑖, are calculated with the limiter function superbee and 𝑔^1𝑖, 𝑔^5𝑖, with the limiter function minmod.

For the evaluation of 𝑔^𝑚𝑖 and 𝑔^𝑚𝑖+1 in (2.4), it is necessary to calculate the spectral decompositions of the conservative variables increments at the interfaces 𝑖−1/2, 𝑖+1/2, and 𝑖+3/2. In the context of three-dimensional not structured meshes of tetrahedrons, the identification of the cells 𝑖 and 𝑖+1 is intuitive (they are two cells that share a face) but the determination of the points 𝑖−1 and 𝑖+2 is not direct. If two tetrahedrons that share a face are analyzed, the nodes not belonging to the common face can be used as representative points for 𝑖−1 and 𝑖+2. Then, these points can be used as imaginary cells. In this work these ideas have been implemented, being the nodal values calculated as a pondered average of the conservative variables between all cells that are in contact with the nodes 𝑖−1 and 𝑖+2. Such pondered average is given by

𝐔_node𝑘=∑^𝑛𝑖=1𝐔_cell𝑖/𝑑_{𝐺𝐶cell𝑖−node𝑘}∑^𝑛𝑖=11/𝑑_{𝐺𝐶cell𝑖−node𝑘},(2.12) where 𝑑_{𝐺𝐶cell𝑖−node𝑘} is the distance that separates the gravity center of the cell i from the node k, and n is the cell number in contact with the node k.

The treatment of the boundary conditions is carried out through the imaginary cells technique [2–4]. Five different types of boundaries are considered (1) subsonic inlet; (2) supersonic inlet, (3) subsonic exit, (4) supersonic exit, (5) nonpenetration (solid boundary and symmetry).

2.2. Reentry Equations of Motion

The choice of a suitable set of coordinates and parameters of the trajectory to describe the movement of an object in atmospheric reentry is inherent to any investigation of guided spacecraft. To analyze a reentry trajectory it is appropriated to describe the motion of the center of mass using a set of elements known as Flight Coordinates [12].

Thus, the flight coordinates are described by the six orbital elements: magnitude of the position vector, r, longitude, 𝜃, latitude, 𝜑, magnitude of the velocity vector, 𝑣, flight-path angle, 𝛾, and heading angle, 𝜓 (azimuth of the velocity). At every moment this object is under the influence of a total force, 𝐹⃗, composed by the gravitational force, 𝐹_𝐺, the aerodynamic force, 𝐴⃗, and the force of propulsion, 𝑇⃗:

𝐹⃗=𝑇⃗+𝐴⃗+𝐹_𝐺.(2.13)

The gravitational force is always present. For nonpowered flight, the propulsion force is zero, while for flights outside the atmosphere, the aerodynamic force vanishes.

To derive the equations of motion, we must use an Earth-fixed reference system. The kinematics equations of motion are [12]

𝑑𝑟𝑑𝜃𝑑𝑡=𝑣sin𝛾,𝑑𝑑𝜙𝑡=𝑣co𝑟sc𝛾osc𝜙os𝜓𝑑𝑡=𝑣cos𝛾𝑟sin𝜓.,(2.14)

It is desirable to separate the aerodynamic force into two components and define the tangential component of the nongravitational force, 𝐹_𝑇, along the velocity vector, and the normal component, 𝐹_𝑁, orthogonal to the velocity at the aerodynamic plane. When we have in plane flights, the normal vector 𝐹_𝑁 is in the plane (𝑟,⃗𝑣)⃗, the vertical plane, and there is no lateral force. However, it is possible to create a lateral component of this force, which has the effect of changing the orbital plane. The non-gravitational force is then decomposed into a component on the vertical plane and orthogonal to the velocity vector, and a component orthogonal to this plane using the bank angle, 𝜎.

The force equations are [12]:

𝑑𝑣𝑑𝑡=𝑚1𝐹^𝑇−1𝑚𝐹^𝐺𝑣𝑑𝛾sin𝛾+𝜔²𝑟cos𝜙(sin𝛾cos𝜙−cos𝛾sin𝜓sin𝜙),𝑑𝑡=𝑚1𝐹^𝑁cos𝜎−1𝑚𝐹^𝐺cos𝛾+𝑣²+𝜔₂𝑟cos𝛾+2𝜔𝑣cos𝜓cos𝜙𝑣𝑑𝜓𝑟cos𝜙(cos𝛾cos𝜙+sin𝛾sin𝜓sin𝜙),𝑑𝑡=𝑚1𝐹^𝑁sin𝜎−𝜔cos𝛾−𝑣𝑟²cos𝛾cos𝜓tan𝜙+2𝜔𝑣(tan𝛾sin𝜓cos𝜙−sin𝜙)₂𝑟cos𝛾cos𝜓sin𝜙cos𝜙.(2.15) Here 𝜔 is the rotation of the Earth that appears because we have to consider a reference system fixed on the planet.

2.3. Attitude Equations

Knowing the attitude of a space object means knowing the orientation of an axis system connected to the vehicle related to a vertical reference system. To specify the orientation of a rigid body in space, three independent parameters are needed. These parameters are commonly known as roll, pitch and yaw, the Euler angles.

However, the use of the Euler angles to compute the attitude evolution of a spacecraft is limited: the equations of motion in attitude have singularities for certain values of the pitch angle, namely ±𝜋/2. This limitation was solved with the substitution of the Euler angles with a set of variables known as quaternions [13].

The quaternions are denoted as 𝐪=(𝑞₁,𝑞₂,𝑞₃,𝑞₄). The components of 𝐪 are defined in terms of the Euler angles using the convention 𝑥𝑦𝑧 (or sequence 321) which is found in Goldstein [14].

Kinematics equations of motion have the following form [14]:

𝑞₁̇=12𝜔³𝑞₂̇𝑞₂−𝜔₂𝑞₃+𝜔₁𝑞₄,=12−𝜔³𝑞₃̇𝑞₁+𝜔₁𝑞₃+𝜔₂𝑞₄,=12𝜔²𝑞₄̇𝑞₁−𝜔₁𝑞₂+𝜔₃𝑞₄,=12−𝜔¹𝑞₁−𝜔₂𝑞₂−𝜔₃𝑞₃.(2.16)

The attitude dynamic equations of motion express the temporal dependence of the angular velocity related to the applied torques:

𝐈𝑑𝜔⃗𝑑𝑡=𝑁−𝜔⃗×𝐈𝜔⃗,(2.17) where I is the inertia matrix, and 𝑁 the aerodynamic torque. Thus, with 𝐼₁, 𝐼₂, 𝐼₃ corresponding to the moments of inertia about the main axes of the vehicle, the dynamic equations of attitude are [13]: 𝜔̇₁=𝑁¹+𝐼₂−𝐼₃𝜔₂𝜔₃𝐼₁𝜔̇₂,=𝑁²+𝐼₃−𝐼₁𝜔₃𝜔₁𝐼₂𝜔̇₃,=𝑁³+𝐼₁−𝐼₂𝜔₁𝜔₂𝐼₃.(2.18)

2.4. Aerodynamic Forces

If the vehicle under consideration operates on a symmetry condition, the velocity vector defines this plane of symmetry. So the attitude of the object is properly described by the attack angle, 𝛼, which is the angle between the velocity vector relative to the atmosphere and a vehicle's baseline, normally the longitudinal axis.

The aerodynamic force is decomposed into two components: the force opposite to the direction of motion, called Drag -and part of the non-gravitational force 𝐹_𝑇 in (2.15), and the orthogonal component, called Lift (part of 𝐹_𝑁):

𝐷=12𝐶^𝐷𝐿=1𝜌_atm𝑆_sat𝑉²,2𝐶^𝐿𝜌_atm𝑆_sat𝑉².(2.19) Here, 𝐶_𝐷 and 𝐶_𝐿 are called drag and lift coefficients, respectively, 𝜌_atm is the atmosphere’s density, 𝑆_sat is the satellite’s area and 𝑉 is the relative velocity between the space object and the atmosphere. For each object it is necessary to calculate specific coefficients at every moment of the trajectory. The code developed in the Department of Aeronautics of the UNC [1] evaluates these parameters for the conical object into consideration.

The first phase of the calculation is to identify the object’s surface areas. So, the nodes of the tetrahedrons, whose faces form the surface of the object, are identified in the code input file.

The forces acting on each face have a magnitude equal to the pressure divided by the area. The direction is the incoming normal to the surface, and the point of application corresponds to the geometric center of the face. In this way, the resultant force on the body is calculated as the vector sum of all forces acting on the discretized surface. Finally, the resulting force components in directions parallel and perpendicular to the flow velocity vector are projected, in order to obtain the drag and lift coefficients, respectively.

3. Results

The numerical simulations were performed using as a core calculation code the one developed at UNC [1]. For delineation of the geometric bodies and meshing, we used the application GID 8.0 with temporary license issued by the manufacturer.

GID has been developed as an interface for geometric modelling, meshing, income data, and display results of all types of numerical simulation programs. The different menus can be modified according to specific user needs. The graphical interface adapted to the code is designed to allow the entrance of initial conditions, to mark the object's surface for the calculation of forces, to define the number of iterations, and other parameters inherent to the code. The ultimate objective of the use of GID is writing the data file that enters the UNC's code and the subsequent display of results, from reading the output file written by the code.

The limit of iterations in the code can be determined by the number of steps, or the limit time independently. When it comes to any of these, the program completes its implementation.

The motivation of the simulations performed is to obtain the main aerodynamic characteristics (drag and lift coefficients) of a conical body. They are used in the calculation of the trajectory during reentry into the atmosphere [15]. Due to the fact that the original code did not have the calculation of the pressure forces on a body, the necessary sentences to perform a simple summation of forces on the faces of the tetrahedrons lying on the object were written. For the identification of these faces, the graphical interface developed in GID is used, which has an option to mark the surface of the body.

The angle of the cone (10°) was chosen arbitrarily, in order to guarantee that this value is small. The cone's length (1 m) was chosen to facilitate the calculations.

It is known that the aerodynamic characteristics of a body with the given geometry depend, in the case of nonviscous flow, only on the Mach number and the incidence angle of the free flow, without considering heat transfer or changes in the properties of air. That's the reason why these variables are used as independent ones.

Two sets of simulations were performed. The first case is for a cone-shaped object that enters in the atmosphere by sharp (narrow) side forward, while the second one was calculated in the assumption that the vehicle moves by wide (obtuse) side forward. Both calculations were made for several attack angles and Mach numbers.

It is important to note that the pressure distribution around the body depends only on the free flow Mach number. Then to obtain the 𝐶_𝐷 and 𝐶_𝐿 coefficients, it is necessary to modify only the free flow Mach number. For this reason the velocity is modified by keeping pressure and density constant and equal to those used in the remaining cases. The speed in the z direction is maintained equal to zero so there is flow symmetry around the 𝑥𝑦 plane.

3.1. Sharp Side Forward Case

The dimensions of the volume meshing are length in direction x = 2.2 m, length in direction y = 3 m, distance to surface point of impact in the axis x = 0.2 m, angle of revolution = 180°. The characteristics of the mesh are number of tetrahedrons around 130,000 and number of nodes around 29,000.

The interface with GID requested that a new mesh was drawn every time a new calculation was made, because it needed a new data file for every case considered. For this reason, and because a nonstructured mesh was used, every case had a slightly different number of elements, due to arbitraries processes during meshing.

The criterion used for the meshing is based on

(i)using the symmetry of the flow on a plan to reduce the volume of mesh, and on two planes in the case of zero attack angle,(ii)concentration of all elements near the surface of the cone, in particular areas considered critical for the calculation: impact point and base,(iii)significant reduction of elements in remote areas of the body, to reduce computational cost.

Figures 1(a) and 1(b) display some examples of results obtained for the pressure field over a conical object considering two combinations of Mach number and attack angle.

(a)

(b)

Figure 1 

Pressure distribution on the object: (a) Mach = 1, attack angle = 0°; (b) Mach = 4, a.a. = 15°.

The implemented numerical scheme has the capacity to simulate compressible flows in subsonic, transonic, and supersonic regimens. From Figure 1(a), it is possible to note, for Mach number equal 1, the concentration of sonic waves near of the body edge; however, the shock wave has not yet been formed. The shock wave can be seen clearly in Figure 1(b) when the free Mach flow is 4. Furthermore Figure 1 shows that the pressure distribution, including at the cone base, is different for transonic from that of supersonic flows, this explains the different behaviors of the 𝐶_𝐷 and 𝐶_𝐿 coefficients for different regimes.

Figure 2 shows the drag coefficients for the simulated cases, while Figure 3 shows the lift coefficients. In all cases it is considered the cone angle equal to 10°. The reference area used is the front section of the cone. It is important to note that the viscous effects are neglected to obtain theses figures.

Figure 2 

Sharp Side Forward: Drag coefficient as a function Mach for different angles of attack.

Figure 3 

Sharp Side Forward: Lift coefficient as a function Mach for different angles of attack.

Note from Figure 2 that drag coefficient is increasing accordingly as the attack angle of the cone is increasing too, this phenomenon occurs until 75°. From Figures 2 and 3, it is possible to observe that both coefficients, drag and lift, have lower variations for height Mach numbers and the greater variations occur at subsonic and transonic flows. For an attack angle of 45°, it is produced that the highest lift and the lower occur for an angle of 90°, for this test the lift is negative except for very reduced Mach number.

Figure 4 show some of the results of the reentry simulations. Figure 4(a) plots the variation of attack angle during the reentry trajectory for initial values of this angle corresponding to 0°, 15°, and 30°, while Figure 4(b) presents the trajectory, meaning the altitude as a function of the time, for these angles.

(a)

(b)

Figure 4 

Sharp Side Forward: (a) Attack angles depending on the time; (b) Altitude. Initial attack angles: 0° (blue), 15° (red) and 30° (green).

We can notice the moment when the objects sense the presence of the atmosphere (around 1500 seconds) and begin to experience high variation of the attack angle during their descent. The first trend is to decrease the initial attack angle. However, although the atmosphere promotes an oscillatory movement around the zero attack angle condition, this is very unstable for objects reentering sharp side forward, and it is possible that for certain initial conditions, the atmosphere turns the object and reverses its attitude.

In what concerns the time evolution of the altitude, it could be seen that it is not affected until the height of 150 km. This is a limitation of the atmospheric model that has information on the atmospheric density just until this height. Since the aerodynamic forces depend directly on this parameter, it is assumed that the atmosphere is so faint over 150 km that it does not affect the trajectory. Anyway, under this limit, all trajectories show bouncing movements, indicating that the attitude variation affects the objects dynamics.

3.2. Wide Side Forward Case

In this case the characteristics of the mesh are number of tetrahedrons around 160,000 and number of nodes around 32,000. An additional criterion used for the meshing is based on the concentration of all elements over the cone’s surface, in particular areas considered critical for the calculation: impact point and base. Also smaller elements were used in the zones nearby the cone base and along the axis of the cone in order to determine more precisely detached shock waves. Apart from that, the mesh was made thinner over the sides of the cone compared to the aerodynamic shadow in order to more accurately determine the expansion waves.

To simulate the atmospheric reentry of a cone shaped object by wide side forward, it is takes into account that this object has an attack angle around 180º. Some of the results found for this condition can be seeing on Figures 5 and 6.

Figure 5 

Wide Side Forward: Drag coefficient as a function Mach for different angles of attack.

Figure 6 

Wide Side Forward: Lift coefficient as a function Mach for different angles of attack.

Figure 5 shows the drag coefficients as a function of the Mach number for attack angles 105°, 120°, 135°, 150°, 165°, and 180°. It can be seen that this parameter presents higher values on the condition wide side forward compared with the same attitude angles on a sharp side forward situation for attack angles around the stability condition (180° in this case and 0° in the sharp side forward) as expected. For attack angles equal to or bigger than 135° (corresponding to 45° in the other case), the drag coefficient is very similar. On the other side, the lift coefficients tend to be smaller or maintain the same values (Figure 6).

The result of these differences can be noticed in Figures 7(a) and 7(b). For a conic object that reenters the terrestrial atmosphere wide side forward, small angles of attack (here around 180°) tend to preserve their magnitude during the reentry trajectory as corresponding to an equilibrium configuration. However, when the initial attitude is as far from 180° as 30°, equilibrium is lost during reentry and we found again the oscillations around the stability condition.

(a)

(b)

Figure 7 

Wide Side Forward: (a) Attack angles depending on the time; (b) Altitude. Initial attack angles: 180° (dark blue), 175° (dark red) and 150° (dark green).

The most important consequence of the wide side forward reentry is over the altitude evolution. It can be seen in Figure 7(b) that all the cases present nearly ballistic paths. Although the far from the 180° initial attack angle condition, the less ballistic is the trajectory inside the dense part of the atmosphere, as can be seen by the dark green line that corresponds to an initial attack angle of 150°, the trajectories of initial attack angles equal to 180° and 175° are superimposed.

4. Conclusions

In this paper, a series of numerical simulations have been conducted using a code developed at UNC [1, 2], which solves the Euler equations by the method of finite volume, using unstructured 3D tetrahedral meshes. The code implements a new technique on the introduction of limiting functions, which aims to decrease the numerical viscosity, that is, increasing the contact discontinuities capture accuracy without loss of robustness regarding other TVD methods.

The goal was the calculation of aerodynamic characteristics of a cone under the effects of different flight conditions, to attach the results to other code that calculates the dynamic of atmospheric reentry [7].

It was possible to perform simulations of reentry trajectories for a specific cone under the influence of various aerodynamic effects for a variety of initial attack angles. Simulations show, as expected, that the trajectories are more affected when the object has initially a sharp side forward configuration.

From the obtained results, was compared the numerical slope of the normal force coefficient with analytical results. In none of the cases the errors exceed 6%. This good accuracy of the numerical results permits to induce that the error by not considering viscous effects in the calculation of the aerodynamic coefficients is low, and the obtained trajectories and attack angles evolution during the reentry are reliable.

Although the methodology implemented has been shown to be suitable for calculating the reentry trajectories inside the terrestrial atmosphere, it will be improved with the inclusion of viscous effects in the simulation of the aerodynamic flow in future works.