End-to-end Mars entry, descent, and landing modeling and simulations for Tianwen-1 guidance, navigation, and control system

On May 15, 2021, the Tianwen-1 lander successfully touched down on the surface of Mars. To ensure the success of the landing mission, an end-to-end Mars entry, descent, and landing (EDL) simulator is developed to assess the guidance, navigation, and control (GNC) system performance, and determine the critical operation and lander parameters. The high-fidelity models of the Mars atmosphere, parachute, and lander system that are incorporated into the simulator are described. Using the developed simulator, simulations of the Tianwen-1 lander EDL are performed. The results indicate that the simulator is valid, and the GNC system of the Tianwen-1 lander exhibits excellent performance.


Introduction
The Tianwen-1 mission to Mars was launched on July 23, 2020. On May 15, 2021, after approximately 9 min of the entry, descent, and landing (EDL) phase, the Tianwen-1 lander successfully touched down on the surface of Mars. The Mars EDL phase began at the atmospheric interface at a radial altitude of 125 km and ended with a surface touchdown [1,2]. This paper describes the EDL modeling and simulation of the Tianwen-1 EDL mission.
The Tianwen-1 EDL modeling and simulation are aimed at designing and validating the lander guidance, navigation, and control (GNC) system, which ensures accuracy and safe touchdown. In the EDL phase, the Tianwen-1 lander sequentially performs trim wing deployment, parachute deployment, heatshield ejection, and backshell separation in the highly uncertain Mars environment [3]. Thus, the dynamics of the EDL phase are complicated and involve multiple bodies with significant interactions [4]. To accurately model the dynamics of the EDL phase, we develop an end-to-end Mars EDL simulation tool for the Tianwen-1 GNC system, including high-fidelity models of the Mars atmosphere, parachute, and lander system. This simulation tool can be used to assess the GNC system performance and determine the critical mission and lander parameters.
The National Aeronautics Space Administration (NASA) developed two primary Mars EDL simulation tools, based on NASA LaRC's program to optimize simulated trajectories II [5,6] and NASA JPL's dynamics simulator for entry, descent, and surface landing [7,8]. Based on the two simulators, NASA performed Mars EDL simulations for the Mars Phoenix and Mars Science Laboratory and successfully landed them on the Mars surface [9,10]. Several of our simulator models, including the Mars atmosphere environment and parachute dynamics, have been developed with reference to the foregoing Mars EDL flight experiences [11][12][13], with modifications introduced to the Tianwen-1 lander. To satisfy the requirements of parachute deployment, the trim wing of the Tianwen-1 lander is deployed at a specified velocity to secure a zero trim angle of attack that differs from previous Mars landing missions. Accordingly, the dynamics of the Tianwen-1 lander with the deployed trim wing have been modeled in our huangxyhit@sina.com Nomenclature EDL entry, descent, and landing GNC guidance, navigation, and control IMU inertial measurement unit DEM digital elevation map DOM digital orthophoto map rI position vector vI velocity vector q attitude quaternion ω angular velocity vector CA, CN, CZ axial, normal, and side aerodynamic force coefficients, respectively CAP, CNP, CTP axial and normal aerodynamic force coefficients and pitch moment coefficient, respectively CTx, CTy, CTz roll, pitch, and yaw moment coefficients, respectively C dTy , C dTz pitch and yaw dynamic derivative coefficients, respectively KKI, KAOM, KMEF, K Off inflation factor, area oscillation parameter, Mach efficiency factor, and offloading model parameter, respectively simulator. Moreover, we have developed interpolation methods for accurately simulating the aerodynamics, parachute dynamics, and dynamics processes of trim wing deployment, parachute deployment, heatshield injection, and backshell separation. We have also adopted multi-threading and parallel processing technologies to design our simulator. Based on these methods, our simulator simplifies the calculation of the dynamics and enables the implementation of a Monte Carlo simulation with tens of thousands of runs over a brief period. Compared with the extensive research conducted by the United States on Mars EDL modeling and simulation, Chinese scholars have focused more on the GNC technology associated with the Mars EDL [14,15]. Although certain studies have been focused on parachute dynamics [16][17][18], the frameworks are mainly theoretical and not suitable for engineering applications. As described in this paper, we update the existing Mars atmosphere environment and parachute dynamics model and develop several models to construct a complete Mars EDL simulator for the Tianwen-1 lander. With the above equipment, several simulations for evaluating the capability and performance of the Tianwen-1 lander GNC system are performed.

Tianwen-1 lander configuration
The Tianwen-1 lander consists of heatshield, backshell, and landing platform carrying a rover, as shown in Fig. 1. The backshell contains a trim wing and a parachute system, which are sequentially deployed at scheduled periods. The lander is equipped with navigation sensors, including an inertial measurement unit (IMU), a multibeam radar, a Mars imaging sensor, actuators (including 26 attitude control thrusters), and a variable-thrust main engine. The IMU includes gyros and accelerometers, which measure the three-axis attitude angular velocity and nongravitational acceleration of the lander. The multibeam radar measures the distance to the ground and velocity along each beam direction. The Mars imaging sensors include optical imaging sensors and three-dimensional (3D) laser imaging sensors, which are used to obtain images of the Mars surface and identify safe landing sites while avoiding hazards.

Tianwen-1 EDL sequence phase and operation
The EDL sequence profile of the Tianwen-1 lander, which consists of atmospheric entry, parachute descent, and landing phases, is shown in Fig. 2. The Tianwen-1 EDL mission is initiated upon atmospheric entry, which is regarded to occur at an altitude of approximately 125 km with an entry inertial velocity of 4.8 km/s. When the lander decelerates to 2.8 Mach, the trim wing is deployed, and the angle of attack is trimmed to zero, which is suitable for parachute deployment. When the lander decelerates to 1.8 March, the parachute is ejected from the rear of the backshell, and the parachute descent phase is initiated. In this phase, the heatshield is jettisoned at a specific time after the parachute deployment, and the lander legs are deployed. Once the lander legs are deployed, the landing radar begins to measure the vehicle altitude and velocity with respect to the Mars surface. The parachute descent continues until the navigated altitude and velocity satisfy the requirements of powered descent. Next, the lander releases the backshell and ignites the main engine, and the landing phase begins. In this phase, the lander performs velocity reduction, backshell evasion maneuvers, and hazard avoidance. Finally, it performs soft touchdown on the selected landing site. The Tianwen-1 EDL simulation must model all the EDL phases and operations, which are triggered by the GNC software included in the simulation.

Tianwen-1 EDL end-to-end simulator
The architecture of the simulator is illustrated in Fig. 3

Modeling overview
The six modules described in Section 3 are modeled in the developed Tianwen-1 EDL simulator. This study is focused on the modeling of the Mars environment, multibody dynamics, and GNC algorithm, which are unique and crucial for the Tianwen-1 EDL mission. As described in Section 2, the Mars EDL phase consists of atmospheric entry, parachute descent, and landing phases. In the atmospheric entry phase, the Tianwen-1 lander with the backshell, heatshield, and packed parachute is treated as a rigid body (known as the entry capsule). The entry capsule is subjected to gravity, aerodynamic force, and engine thrust. Modeling the aerodynamic force is a key element of the Tianwen-1 EDL simulator because knowledge regarding the Mars atmospheric environment is insufficient; thus, its atmosphere must be regarded as highly uncertain.
Another key problem is the modeling of parachute dynamics involving the phases of parachute ejection, parachute line stretching, parachute inflation, and stabilized descent. In these phases, for simplicity, the parachute and lander are modeled as two rigid bodies connected by parachute lines. Once the heatshield is ejected, the heatshield dynamics are modeled and solved to verify whether the heatshield is safely separated. The influence of heatshield separation is evaluated based on its impact on the radar beam. Similarly, after backshell separation, the parachute force is exerted on the backshell by parachute lines, and the backshell dynamics are modeled and solved to validate and evaluate the parachute-backshell avoidance. The GNC module provides accurate and reliable state estimates and key operation triggers as well as implements accurate and reliable trajectory and attitude controls. In this context, it is a core module validated by the simulator.

Tianwen-1 lander dynamics and modeling
The lander dynamics are modeled using two reference frames: the Mars-centered J2000 inertial frame and lander body frame. The lander dynamics are described in the Mars-centered J2000 inertial frame as where r I and v I are the lander position and velocity vectors in the Mars-centered J2000 frame, respectively; g mI is the gravity vector in the J2000 inertial frame, which can be calculated according to the lander position and Mars gravity field model contained in the Mars environment module; a thr , a atm , a para , and a hea denote the accelerations caused by the engine thrust, aerodynamic force, parachute drag force, and recoil force of heatshield ejection, respectively; q is the lander attitude quaternion describing the rotation from the J2000 inertial frame to the lander body frame; ω is the angular velocity corresponding to the lander attitude; I is the inertial matrix of the lander; and T thr , T atm , T para , and T hea are the torques associated with the engine thrust, aerodynamic force, parachute drag force, and recoil force of heatshield ejection, respectively. The accelerations a thr , a atm , a para , and a hea are calculated as where m is the lander mass with the consumption ratė m = n i=0 ∥F thri ∥ Ispi ; n is the thrust number; F thri and I spi are the thrust and impulse of the ith thruster, respectively; F atm , F para , and F hea are the aerodynamic force, parachute drag force, and recoil force of heatshield ejection, respectively.
To solve the dynamics described in Eqs. (1) and (2), the forces (i.e., F thri , F atm , F para , and F hea ) and torques (i.e., T thr , T atm , T para , and T hea ) must be calculated. Using the thrust model, the engine force, F thri , can be calculated according to the thruster pulse width commands controlled by the GNC module. Torques T para and T hea are calculated using the parachute dynamics model and heatshield dynamics model, respectively. In particular, F para = 0 and T para = 0 in the atmospheric entry phase and landing phase, and F para and T para have to be solved only in the parachute descent phase for the lander, as described in the following section. The variation of F hea and T hea with time occurs over a period of only 100 ms; they are modeled using the heatshield dynamics module. For the entire EDL phase, F atm and T atm are calculated and modeled as follows: where C q is the rotation matrix from the body frame to the J2000 frame; q p = 1 2 ρv 2 r is the aerodynamic pressure; ρ is the atmospheric density at a current position (provided by the Mars atmosphere model); v r is the airspeed magnitude calculated as the difference in the lander velocity relative to the Mars surface and wind speed; ; C A , C N , and C Z are the axial, normal, and side aerodynamic force coefficients, respectively; C Tx , C Ty , and C Tz are the roll, pitch, and yaw moment coefficients, respectively; and C dTy and C dTz are the pitch and yaw dynamic derivative coefficients, respectively.
The aerodynamic force, moment, and dynamic derivative coefficients can be calculated by interpolating the corresponding aerodynamic coefficient databases [12]. To support the EDL phase, several aerodynamic coefficient databases are modeled according to different phases and configurations based on empirical or simulation data. When the altitude exceeds 44.8 km, the aerodynamic force and moment coefficients (known as static aerodynamic coefficients) are modeled for a specific range of altitudes and total attack angles. When the altitude is less than 44.8 km, and the trim wing is not deployed, the static aerodynamic coefficient is modeled for different Mach numbers and total attack angles. In the trim wing deployment process, the static aerodynamic coefficient is modeled according to the deployed angles of the trim wing and total attack angles. After trim wing deployment, the static aerodynamic coefficient is modeled considering the Mach numbers, attack angles, and sideslip angles. In particular, during the heatshield and backshell separation process, the separating object moves extremely close to the lander, and the two entities exert considerable aerodynamic effects on each other. Therefore, the static aerodynamic coefficient is modeled considering the Mach numbers, attack angles, sideslip angles, and distances between the lander and separating object. The aerodynamic derivative coefficients are modeled considering the Mach numbers, attack angles, and sideslip angles, accounting for the cases of trim wing deployment and nondeployment.
Using these aerodynamic coefficient databases, the orbit and attitude dynamics can be simulated throughout the Mars EDL phase, especially operations. Because of the uncertainty in aerodynamics, aerodynamic coefficient bias databases are constructed in the simulator. These are inconsistent with the corresponding coefficient databases and are used in the Monte Carlo analysis, as described in the next section.

Atmosphere and wind modeling
Atmospheric and wind modeling is performed to specify the atmospheric temperature, density, and wind velocity calculated by interpolating the Mars atmospheric environment databases according to the altitude in the simulator. The databases contain nominal values and their bias ranges. The nominal values and bounds of density, temperature, and horizontal wind speed are shown in Figs. 4-6, constructed with reference to the Mars Science Laboratory (MSL). The vertical wind speed is assumed to be zero when the altitude exceeds 20 km and set as a zero mean value with a 20 m/s deviation when the altitude is less than 20 km [9]. These bounds represent the uncertainties in the atmosphere, which are considered in the Monte Carlo analysis described in the next section.
After determining the atmospheric temperature, density, and wind velocity, the Mach number, aerodynamic pressure, and angles of attack and sideslip can be calculated. These values are inputted to the orbit and attitude dynamics model to solve the aerodynamic force and torque.

Parachute dynamics and modeling
The Tianwen-1 lander parachute is similar to that of the MSL described in Ref. [13] and contains a single riser and four bridles. The parachute is connected to the backshell by the four bridles, and the parachute drag force is calculated according to the change in the length of each bridle.
Parachute descent consists of the following phases: parachute ejection, parachute line stretching, parachute inflation, and stabilized descent. Before the parachute lines stretch, solving the parachute dynamics is unnecessary, and only the parachute ejection force must be considered because the recoil force of ejection is imposed on the lander. The parachute ejection phase lasts from mortar firing at time t 0 to the parachute bag reaching the mortar muzzle velocity at time t open . In the parachute ejection phase, the mortar ejection force is set as 100 kN for Tianwen-1. Assuming that the parachute lines stretch at time t ls , and backshell separation occurs at time t sepB , the parachute drag force, F para , on the lander (Eq. (2)) can be calculated. This force is estimated according to the change in the length of each bridle, which can be determined from the poses of the parachute and backshell. After calculating F para , the corresponding torque, T para , on the lander can be calculated according to the position of the connecting point between the bridles and backshell.
After the parachute lines stretch, the parachute dynamics are modeled similar to that of the Tianwen-1 lander. In the model, the forces on the parachute are gravity, aerodynamic force (F atmP ), and parachute drag recoil force (F paraP ). The recoil force is calculated as and the corresponding torque, T paraP , on the parachute can be calculated according to the position of the connecting point between the bridles and riser.
The aerodynamic force (F atmP ) and torque (T atmP ) acting on the parachute are modeled as where C qP is the rotation matrix from the parachute body frame to the J2000 frame; q pP is the dynamic pressure acting on the parachute; D refP is the reference diameter with S refP = πD 2 refP 4 ; K KI , K AOM , K MEF , and K Off denote the inflation factor, area oscillation parameter, Mach efficiency factor, and offloading model parameter, respectively (details are found in Ref. [13]); and C AP , C NP , and C TP are the axial and normal aerodynamic force coefficients and pitch moment coefficient, respectively, describing the static aerodynamic coefficients of the parachute because the parachute is axisymmetric.
The static aerodynamic coefficients of the parachute, C AP , C NP , and C TP , can be calculated by interpolating the corresponding parachute aerodynamic coefficient databases in our simulator; they are modeled for the total angle of attack, as shown in Figs. 7-9. The parachute dynamics solution step must not exceed 1 ms due to the considerable dynamic oscillation of the parachute.

Heatshield and backshell dynamics modeling
The heatshield or backshell dynamics are solved when the objects separate from the landing platform, and the model of the orbit and attitude dynamics is similar to that of the lander described in Eq. (1). In the dynamics model, the forces acting on the heatshield are gravity and aerodynamic force; the forces acting on the backshell are gravity, aerodynamic force, and parachute drag force. The aerodynamic force is calculated by interpolating the corresponding aerodynamic coefficients of the simulator, and the parachute drag force is calculated according to the change in the length of each bridle, as described in Eq. (5). However, the heatshield is ejected through the explosion of four bolts, and the explosive force acting on each bolt is imposed on both heatshield and landing platform. The total explosive force impact on the heatshield is defined as F heaH = −F hea = 4 i=0 F heaHi over a period of only 100 ms. The explosive force direction is along the x axis of the lander, and its magnitude is modeled as a function of the distance between the heatshield and landing platform, as shown in Fig. 10. The torque generated by the explosive force can be calculated according to the installation location of the explosive bolts.

GNC module for Tianwen-1 EDL
The GNC module contains eight GNC modes: attack angle trim mode, lift control mode, parachute descent control mode, powered deceleration mode, hover and imaging mode, hazard avoidance maneuver mode, slow descent mode, and no-control mode. In each mode, the GNC module processes data from sensors, performs real-time GNC calculations, and transmits the control signals for orbit and attitude controls. The details of the Tianwen-1 lander GNC algorithm are presented in Ref. [3].

Mars EDL simulation results under nominal conditions
This section describes the simulation under nominal conditions. The initial states of the lander at the atmosphere interface are set as those used in the actual flight operation on May 15, 2021, as listed in Table 1.
The initial position and velocity are expressed in the Mars-centered J2000 frame.  The angular velocity is less than 0.1 • , which satisfies the requirements of safe and soft landing. In this simulation, note that the duration of the parachute descent phase is between 273 and 420 s. During this period, parachute inflation and area oscillation can cause high angular rate oscillation. Meanwhile, rate-damping controllers are used for the pitch and yaw channels with long deadbands, and the angle and angle rate error are controlled for the roll channel. As a result, the pitch and yaw angle rates oscillate during the parachute descent phase, as shown in Fig. 14.
The Mach number, angles of attack and sideslip, and

Monte Carlo analysis
To assess the performance and robustness of the GNC scheme designed for the Tianwen-1 lander, a Monte Carlo analysis is performed using the constructed simulator. The analysis is implemented considering the uncertainties in the key parameters and coefficients of the environment and lander as well as the uncertainties in the initial states of the lander. The uncertainties in the atmospheric parameters and aerodynamic coefficients are described in Section 4, and the uncertainties of the initial states and characteristics of the lander are listed in Table 2 (specified as 3σ or min/max values dispersed from the nominal values). The uncertainties in position and velocity are expressed in the orbit coordinate system, in which the xy plane coincides with the orbit plane, the x axis is along the radial direction, the y axis points to the velocity direction, and the z axis is along the normal direction of the orbit plane. The uncertainties in the sensors and actuators are also included in the simulation; however, they are not presented in this paper.    pressure is within 560 ± 200 Pa, and the altitude is between 4 and 15 km. At parachute deployment, the vertical velocity lies within −1.5 ± 0.6 m/s, the horizontal velocity is less than 0.9 m/s, the angular velocity is less than 1 ( • )/s, and the angle between the x axis and vertical direction is less than 2.5 • at touchdown. The Monte Carlo simulation results indicate that the constraints are satisfied in all 10,000 iterations.
In addition, the distributions of the propellant consumption and duration of EDL are examined to analyze the landing mission safety. The results shown in   The Monte Carlo simulation generates a considerable amount of lander information during EDL that can be used to fully assess the performance of the GNC system. Not all information is presented here because of limited paper length.

Conclusions
An end-to-end Mars EDL simulator for the Tianwen-1 GNC system is developed. The multibody dynamic models for the lander, parachute, backshell, and heatshield in the simulator are described. Using the developed simulator, simulation under nominal conditions and Monte Carlo statistical analysis are performed to assess the performance and robustness of the Tianwen-1 lander GNC system. The simulation results indicate that the simulator is valid, the Tianwen-1 lander GNC system for the Mars EDL exhibits high performance, and the key parameters in the EDL phase satisfy the requirements. The simulator plays a key role in the development of the Tianwen-1 lander GNC system and can be applied to the next generation of Mars exploration systems.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.