Use of Modern Digital Software to Model the Motion Dynamics of Spacecraft during Landing

Landing systems for future space missions in Earth and Mars require trustable technologies capable of achieving their aim in the most accurate way possible. For this purpose, systems should go through rigorous dynamic simulations run by precise and efficient software. This study aims to approximately determine the dynamic motion of a landing vehicle using the modern digital software of Universal Mechanism and MATLAB. Universal Mechanism applies the classic mechanics theory on a model based on the geometry of the spacecraft, taking into account the environmental conditions that affect its motion and the properties of the ground to resist its impact [1]. The forces implied in the vehicle phase of descent are also included in MATLAB code to calculate the landing area of the vehicle according to its re-entry velocity [2-4]. Studies were conducted for different initial conditions and approaches to the surface. As a result, the values of the arising overloads and forces acting on the descent vehicle were obtained [1]. The data provided by the simulations conclude the safest landing options that should be taken into account for the success of future missions.


Introduction
Modern digital software is essential in every engineering field nowadays. Both public and private companies require some kind of guarantee that their technical projects are going to be successful, since funds invested in them are too high to take any risk.
The space sector is undoubtedly the most expensive sector of all. The slightest mistake or miscalculation can be disastrous, leading to a loss of over 1.3 billion of euros (in the case of the ExoMars missions [5]) and very valuable time. To avoid this situation, every system is thoroughly tested in both digital and physical simulations under all conditions and eventualities.
For instance, in the United States, landing systems are designed, tested and evaluated using end-to-end high-fidelity computer flight simulations and end-to-end real-time hardware-in-the-loop simulation test bed, based on NASA Langley Research Center`s Program to Optimize Simulated Trajectories II and NASA Jet Propulsion Laboratory's Dynamics Simulator of Entry, Descent and Surface Landing (DSENDS), respectively [6].
Generally, scientifically interesting landing sites are not flat and contain many landing hazards including significant variation in terrain elevation, craters and rocks [7]. The touchdown spot of previous missions could only be pre-specified to lie within an uncertainty ellipse with a semi-major axis of several kilometers, ranging from 200 to 300 km in the cases of the Pathfinder and Viking landers, and to 80 km in the case of Mars Exploration Rovers. Therefore, future landing systems must possess a higher degree of autonomy than previously used in order to carry out a successful landing at these difficult to reach sides [7][8].
To improve the accuracy of the landing, new navigation systems have been developed over the years. Such is the case of Magnolia-1, a robust algorithm for autonomous hazard detection and avoidance (HDA) that requires no input from classical inertial sensors, since the hazard assessment and navigation is done by using imaging camera [9]. Hypersonic guidance has also been studied as an option. However, the Small Divert concept (ballistic entry followed by subsonic parachute deployment) does not require it or real-time terrain recognition. Instead, it relies on a priori orbital observation to identify small, multiple safe zones within a larger hazardous region and additional descent propellant to land at the fueloptimal safe zone [7]. Conversely, the Autonomous Landing Hazard Avoidance Technology (ALHAT) system module is a highly advanced integrated sensor suite that enables navigating safely to a landing site with a diameter of 15m and within 100m of a previously specified location [8,10].
This study focuses on the virtual simulation of a landing vehicle in its approach to the surface and during its touchdown on Earth and Mars when there is no parachute deployment. The dynamic motion of the vehicle is estimated and analyzed using modern digital software Universal Mechanism and MATLAB. With the resulting data, vehicle and mission designers will be able to select the optimal landing parameter sequence provided according to the characteristics of the vehicle and the surroundings.

Descent Vehicle Model Used for the Simulations in Universal Mechanism
The descent vehicle can be constructively represented as a system of two completely rigid bodies: the device itself and the payload inside of it. The payload is connected to the vehicle by a passive damping system, which means that it is not isolated and secured for the vibrations that may receive from the main body when the impact takes place. The parachute system compartment has been withdrawn, leaving the geometry that is presented in Fig. 1. This change in the shape of the device only affects the length measure, which now has a value of 1.597 m. The rest of the measures are the same as in the model used for calculations in MATLAB (see Fig. 2).
In terms of its inertial characteristics, the mass of the device is 964.9 kg and its moments of inertia in the system of coordinates attached to the descent vehicle have a value of 0.303 kg•m2, 0.237 kg•m2, 0.252 kg•m2 respectively.
With regard to the payload, it has been depicted as a container with a non-holey toroid surface. Its mass has a value of 229.6 kg and its moments of inertia in the system of coordinates attached to the descent vehicle are 0.025 kg•m2, 0.033 kg•m2, 0.027 kg•m2 respectively.

Landing Simulations Using the Software Universal Mechanism
Equations Universal Mechanism simulates kinematics and dynamics of planar and spatial mechanical systems, representing them as systems of rigid bodies connected by various kinematical pairs and force elements [11].
The landing is simulated in Earth (g = 9.81 m/s2) and Mars (g = 3.73 m/s2) for two different types of soil: boulder and sand. Their properties to resist the impact of the vehicle and friction are modeled with the four coefficients: A, B, C and μ. To find the first three coefficients, it is necessary to have available data about the resistivity for various velocities and its approximation by a square trinomial. The coefficients can then be calculated using the data from the penetration depth obtained from the impact of the vehicle for different velocities. Additionally, the friction coefficient (μ) is determined experimentally, since it depends on the nature of the material and the relative speed between the contact surfaces (among other factors).
The simulations are done for different body and ground inclinations, varying them in steps of 5°. The symmetrical axis of the landing vehicle can be deviated from the vertical axis of the planet from 0° to 20°, while the ground can vary its inclination from -20° to 20°. All in all, this study analyses a total of 180 landing simulations.

Landing Calculations Using the Software MATLAB
The landing area of the vehicle is calculated from the dispersion of intersection points between the trajectory of the vehicle and the surface of the planet using MATLAB. MATLAB was chosen for this task because it analyses data iteratively and designs processes with a programing language that expresses matrix and array mathematics directly.
Moreover, it has professional toolboxes that have been rigorously tested to facilitate the efficient development of the code [12]. The object referenceEllipsoid included in the code represents the World Geodetic System of 1984 (WGS 84) reference ellipsoid [13] for calculations on Earth, and the Mars reference ellipsoid for calculations on Mars. Likewise, the function gravityzonal implements a mathematical representation of zonal harmonic planetary gravity based on planetary gravitational potential (excluding the potential due to planet rotation) [14] for both Earth and Mars models.
Lift and drag depend on the shape of the landing vehicle, the angle of attack, the velocity and the atmospheric characteristics of the planet. The graphs that show the variation of the lift and drag coefficients depending on the angle of attack and the velocity of the body [2,15] for Mars were digitalized using the online software WebPlotDigitizer (Version 4.2) [16] and the data obtained were implemented in the code. For Earth calculations, however, there was no available data about the lift and drag coefficients. Therefore, landing calculations on this planet took into account an approximated constant value of 0.5 for the drag coefficient due to the conical shape of the descent vehicle [17] and a value of 0 for the lift coefficient.
With regards to the atmospheric parameters, they are calculated (in the case of Earth) using the function atmoscoesa, a mathematical representation of the 1976 COESA United States standard lower atmospheric values [18], and interpolating the experimentally obtained values of Viking 1 profile data provided by NASA [19] (in the case of Mars).
The code starts establishing the initial parameters of the system; such as the initial position, the pitch and the module of the initial velocity of the vehicle (reentry velocity direction was set to be parallel to the Equator), all of them expressed in the geographical system of coordinates. After calculating all the forces implied in the vehicle phase of descent, the differential equation system is solved using the function ode45, which utilizes a medium order method [20] to iterate the code until the final position is achieved (the vehicle intersects the surface of the planet).

Dynamic Motion of the Descent Vehicle Modeled in Universal Mechanism
The results obtained from the Universal Mechanism simulations consist of a set of images that represent the motion of the vehicle and the values of the arising overloads and forces acting on it and on the payload. For each planet, type of ground and inclination of the main body the value of the overload versus time is represented, showing its variation depending on the angle of the landing surface.
What is important to look at when comparing these graphics is the maximum value of the overload and the time it takes to reach zero [1].
Simulations on Earth. Landing on Soft Soil. In Fig. 2, it can be noticed that the maximum value of the overload is found in the 20° surface inclination simulation. However, the results from the simulations with surface inclination angles of -20° and -15° should not be dismissed, since they present a second peak. Therefore, neither of these simulations would represent a safe landing option for the vehicle. Paying attention to the time the overload takes to dissipate, it could be said that the best landing options would happen for surface inclination angles of 5°, 0° and 10°.  Nevertheless, the best landing options should be chosen taking into account the evolution of the overload not only in the vehicle, but in the payload as well. Fig. 3 shows some changes in the values of the overloads compared to Fig. 2. The overload with the greatest amplitude happens for a surface inclination angle of 15°, followed by the results for angles 10° and 20°. The second peaks happen in the same cases as in Fig. 3, so they should be avoided in order to protect the payload from damage.
All in all, for soft soil landings on Earth with no body axis inclination, it could be concluded that the safest landings (for both the vehicle and the payload) would occur for surface inclination angles of 5° and 0°, while the most dangerous ones would occur for angles of -20°, -15°, 20° and 15°.
Simulations on Earth. Landing on Hard Soil. In Fig. 4, it can be clearly seen that the maximum value of the overload is achieved for the 20° surface inclination simulation, since there is a notorious difference between the amplitude of this overload and the amplitude of the rest of overloads from other simulations. However, the results from the simulations with surface inclination angles of -20°, -15°, -10°, -5°, 15° and 0° should not be dismissed, since they present a second peak. Therefore, neither of these simulations would represent a safe landing option for the vehicle. Paying attention to the time the overload takes to dissipate, it could be stated that the best landing options would happen for surface inclination angles of 5° and 10° (besides, they are the only two inclination angles that do not appear to have a second peak in their overload values). Regarding the overload of the payload (see Fig. 5), the second peaks happen in the same cases as in Fig. 5, so they should be avoided in order to protect the payload from damage. The overload with the greatest amplitude happens for a surface inclination angle of 15°, followed by the results of angles 20° and 10°. Conversely, the fastest dissipation of the overload happens for angles 5° and 15°. Therefore, for hard soil landings on Earth with no body axis inclination, it could be concluded that the safest landings (for both the vehicle and the payload) would occur for a surface inclination angle of 5°, while the most dangerous ones would occur for angles of -20°, -15°, 20° and 15°.
Simulations on Mars. Landing on Soft Soil. For landings on soft soil on Mars, the maximum value of the overload acting on the main body is achieved for the 15° surface inclination simulation (see Fig. 6). However, it should be noticed that the -20° inclination curve slightly deviates from the common shape that the other inclinations have, showing the foundation of a possible second peak. Therefore, these simulations should be discarded as safe landing options. The simulations that take the littlest time to dissipate the overload are the 5° and 0° simulations. The maximum value of the overload acting on the payload, however, is obtained for the 10° surface inclination simulation (see Fig. 7), followed by the one obtained for the 5° surface inclination simulation. In this case, the -20° inclination simulation also shows the foundation of a possible second peak, and the fastest simulations that reach an overload value of zero are also the same. As a conclusion, for soft soil landings on Mars with no body axis inclination, the safest landings (for both the vehicle and the payload) would happen for the surface inclination angles of 0° and 5°, while the most dangerous ones would occur for angles of 15°, 10° and -20°.
Simulations on Mars. Landing on Hard Soil. Fig. 8 reveals that the simulation with the 20° surface inclination reaches the maximum overload value acting on the main body by far, almost tripling the value reached in the soft soil simulation (see Fig. 8). However, the simulations with the -20°, -15°, -10° and -5° angles should also be taken into account as dangerous landings due to the existence of second peaks in their curves. On the contrary, the safest landings would be achieved for the surface inclinations of 5° and 10°, since the overload tends to zero faster in the curves of their simulations. The overload acting on the payload is shown in Fig. 9. The maximum overload value can be seen in the 15° surface inclination simulation, followed closely by the value of the 10° simulation. As it happened for the overload acting on the main body (see Fig. 9), the simulations -20°, -15°, -10°, -5° and 20° have a second peak in their curves that could be very harmful for the payload. Talking about the time it takes for the overload to fade, the best results appear to come from the 5° and 0° simulations. In summary, the safest landing that could be chosen for hard soil on Mars with no body axis inclination would be the landing on the 5° inclined surface. By contrast, the more dangerous landings would happen for -20°, -15°, -10°, -5° and 20° surface inclinations.
It would be more convenient to choose the landing option that adapts more suitable to the descent vehicle structural requirements than to the payload structural requirements, since the overload acting on the payload could be reduced implementing the appropriate damping systems.
After analyzing the presented data, it can be concluded that the surface inclination and type of soil have a notorious impact on the landing, as well as the deviation of the vehicle vertical axis. The greater the difference between the axis deviation and the surface inclination angles is, the more dangerous the landing becomes. For instance, landing with an axis deviation of 5° on a 5° inclined surface provides a less maximum overload value than landing on a 0° inclined surface. Furthermore, it is more likely for second and third peaks to appear in the overload curves of surface inclination angles that differ significantly from the vehicle axis inclination angle.
Landings on hard soil develop greater values of overloads than landings on soft soil, which affect the final conditions of the payload. The landing impact is more catastrophic on hard soil surfaces, since a greater percentage of the impact energy is absorbed by the vehicle. Even though the percentage of the total overload transmitted to the payload is less in hard soil landings, the overload acting on the payload is still greater in this kind of landings than in soft soil landings.
Therefore, the optimal approach options that provide the safest landings would be those that occur on a slightly inclined surface made of soft soil for axis deviations as small as possible.
The value of the overload is greater for landings on Mars than for landings on Earth. Furthermore, the third and second peaks of the curves are less sharp in Mars landing simulations, meaning that the overload value at those points is prolonged in time.
Not all re-entry vehicles can withstand the overload values that might act on them while landing under certain conditions. Some vehicles may be prepared for high maximum overload values but not for medium repetitive maximum overload values or vice versa. In either of those cases, given a certain landing surface, the vehicle axis inclination angle could be adapted to perform the most suitable landing according to the vehicle structural characteristics.
To ensure the success of the mission, designers should have all these considerations into account when proceeding to choose either the best vehicle axis deviation angle or the most adequate materials and structures for the descent vehicle and its corresponding payload.

Dynamic Motion of the Descent Vehicle Modeled in MATLAB
The code created in MATLAB starts establishing the initial parameters of the system; such as the initial position, the pitch and the module of the initial velocity of the vehicle (re-entry velocity direction was set to be parallel to the Equator), all of them expressed in the geographical system of coordinates. After calculating all the forces implied in the vehicle phase of descent, the differential equation system is solved using the function ode45.
The calculations conducted estimate the landing trajectory of the descent vehicle and represents it along with the altitude versus time. The planet where the hypothetical landing takes place, as well as the initial velocity and pitch angle of the vehicle, are varied in order to analyze the influence of these factors in the vehicle motion.
Landing Trajectory and Altitude. Landing on Earth. The trajectory and altitude versus time of the vehicle represented in Fig. 10 and Fig. 11 were calculated for an initial velocity of 4000 m/s and a pitch of -20°. The vehicle, as it can be seen in Fig. 10, follows a parabolic trajectory during its descent motion. At the end of the landing, the trajectory forms an angle of 90° with the landing surface (taking into account the hypothesis settled for the calculations in section 3.2). Landing Trajectory and Altitude. Landing on Mars. Concerning the trajectory and altitude versus time of the vehicle on Mars, they were also calculated for an initial velocity of 4000 m/s and a pitch of -20° (see Fig. 12 and Fig. 13). However, in this case, the trajectory and altitude versus time are pictured with a linear shape due to Martian atmospheric characteristics, since its density varies almost linearly with the altitude.  Furthermore, the density of the Martian atmosphere is around 1% of the density of the Earth atmosphere [21], which means that drag and lift forces barely act on the vehicle compared to the gravity force.
Landing Site Location and Area. Landing on Earth. The estimated landing positions of the vehicle on Earth depending on the pitch angle for an initial velocity of 3900 m/s are located in an area with a radius of 74300 m (see Fig. 14). Conversely, for initial velocities of 4000 m/s and 4100 m/s, the radius of the landing areas measure, respectively, 75300 m/s and 76200 m/s (see Fig. 15 and Fig. 16). Therefore, the possible landing area of the vehicle increases with the initial velocity.   Nevertheless, the greatest difference between landing positions do not happen varying the initial velocity, but varying the re-entry pitch angle. In fact, as it can be seen in Fig. 14,  Fig. 15 and Fig. 16, the distance between the landing positions increases exponentially as the pitch angle approximates to zero.
Landing Site Location and Area. Landing on Mars. Likewise, the results obtained for the landings on Mars show that for an initial velocity of 3900 m/s an area with a radius of 120000 m is originated, enclosing the possible landing positions for the vehicle depending on the pitch angle (see Fig. 17). For higher velocities, the radius of the area increases, leading to a value of 121300 m and 122400 m for the 4000 m/s and 4100 m/s velocities respectively (see Fig. 18 and Fig. 19). Again, the possible landing area of the vehicle increases with the initial velocity.

Conclusion
The dynamic motion of a descent vehicle during its landing stage requires a deep analysis in all the possible situations in which it could be modified to ensure the reliability of the landing systems applied.
This project has studied the influence that the landing surface inclination and type of soil have on the overload acting on the descent vehicle, as well as the influence of the vertical axis deviation of the vehicle in its approach to the ground. The results have concluded that the greater the difference between the axis deviation and the surface inclination angle is, the more dangerous the landing becomes and the more likely it is for second and third peaks to appear in the overload curves. Furthermore, the landing impact is more catastrophic on hard soil surfaces, since they develop greater values of overloads than landings on soft soil.
Designers could avoid these situations by implementing a control system able to either adapt the vehicle vertical axis to a more convenient inclination given a certain landing surface, or to change the initial conditions of the vehicle re-entry trajectory to make it land on a more suitable surface. If the second option is chosen, it would be appropriate to underline that changing the value of the initial velocity does not affect as much the trajectory of the vehicle as it does changing the value of the pitch angle. In fact, the distance between the possible landing positions increase exponentially as the pitch angle approximates to zero.
All in all, the success of the landing mission strongly relies on the ability of the vehicle and mission designers to select the optimal landing parameter sequence provided according to the characteristics of the vehicle and the surroundings.