Yng DONG, Jinzhong DING,*, Chunjie WANG, Hui WANG, Xueo LIU

a School of Mechanical Engineering and Automation, Beihang University, Beijing 100191, China

b State Key Laboratory of Virtual Reality Technology and Systems, Beihang University, Beijing 100191, China

KEYWORDS Mars lander;Attitude control;Uncertain terrain;Stability analysis;Monte Carlo simulation

Abstract Safe soft landing of the lander is vital to the Mars surface exploration mission.Analysis and verification of the landing stability under uncertain terrain play an important role in lander design.However,the effect of uncertain terrain is ignored in most existing studies,making the analysis incomprehensive and increasing the risk of landing failure in practice.In this paper,a Mars lander with 10 attitude control thrusters is introduced and its dynamics model is then established considering plastic deformation parts and nonlinear contact forces. The effectiveness and accuracy of the dynamics modeling method are verified by experiments with an average relative error of 10%.In order to carry out the dynamics simulation with high-fidelity terrain,a terrain generation method based on statistical data is proposed. Through Monte Carlo simulation under a 50 m × 50 m randomly generated landing terrain,the stability of the lander and the effects of attitude control thrusters are analyzed.The results show that the failure rate is 5.5%,in which the primary failure forms are overturning and abnormal attitude. When the landing simulations are repeated without thrusters,the stable ratio decreases from 94.5%to 90.7%,suggesting the positive effects of attitude control thrusters in improving landing stability under rough terrain.

1. Introduction

With the development of deep space technology, one detector after another was launched to extraterrestrial planets,satellites, and small celestial bodies. Soft landing on celestial bodies is important and challenging in deep space exploration missions, especially in missions that need to land rovers safely on celestial bodies with high gravity, such as the Moon and Mars1–4. Recently, the successful landing of Perseverance5and the launch of the Chinese first Mars lander Tianwen-16have attracted extensive attention of researchers and made Mars exploration a research hotspot. Ways to alleviate the vibration and reduce the impact effectively during the landing process are key concerns in designing a lander. Among all the landing schemes, landing with legged landing gear is well-appraised and widely used due to its simple structure and high reliability7–10.

In order to design a high-stability soft landing lander, it is necessary to fully evaluate and optimize its soft landing performance before launching it onto the outer space planet11–12.Physical experiments and dynamics simulations are currently the two most effective means of analyzing the soft landing performance. Although experiments with prototypes have high credibility13–14, the difficulty in non-earth gravity simulation and high prototype manufacturing cost make the experiments only capable of evaluating limited landing cases. On the contrary, the dynamics simulation is cheap, and it has the advantage of simulating complex landing conditions. Therefore, the finite element(FE)method15and the multi-rigid-body method are widely used for soft landing simulations.The FE method is suitable for analyzing highly nonlinear systems,including large deformation, nonlinear material, and contact16. Some researchers established the FE model of landers to analyze the structural response under harsh loading conditions17. The multi-rigid-body dynamics modeling and analysis method assume that the components are rigid and focus on the macro information such as velocity and displacement during landing.Although the information of deformation and stress is ignored,this method has the advantages of fast computation and low computer memory consumption, which make it widely used in configuration optimization,limit condition traversal,stability boundary analysis,and Monte Carlo simulation18–19.However, errors that arisen from simplification should be considered to give a reliable prediction.

In the process of lander modeling,two issues require special attention.One issue is the simulation of contact forces between the footpads and landing terrain, which may significantly affect simulation accuracy. Nonlinear spring damping model20, plastic contact model21, and model based on force–displacement curve22were developed for contact modeling in recent studies. The other issue is the modeling of landing terrain. Generally, researchers simulated the soft landing using an ideal flat plane to simulate the landing terrain. However,the assumption of landing on an ideal flat plane neglecting the influence of rocks and meteorite craters leads to deviations in the results between simulation and landing in practice.Establishing complex landing terrain based on measured data is helpful for comprehensive analysis of landing performance,such as the studies of Rosetta Lander (comet exploration)23–24. However, it is impossible to analyze the landing terrain completely. For landing on Mars, landing stability in rugged terrain can only be guaranteed by landing gear design and automatic control of attitude thrusters.Therefore,the stability of Mars lander under uncertain terrain needs to be studied deeply.

The Mars lander considered in this paper uses the shutdown-at-touchdown25–26landing mode to reduce the landing velocity and impact force. In addition, 10 attitude control thrusters belonging to Guidance Navigation Control (GNC)system are considered.We aim to activate them during landing to improve the terrain adaptability stability. The main contribution of this paper is that an advanced terrain simulation method is proposed for landing performance evaluation. A high-fidelity dynamics model of the legged Mars lander is established, and its accuracy is verified by prototype experiments. Using the proposed terrain generation methods, the effects of uncertainties of the landing terrain are investigated for the first time,which helps estimate the landing performance comprehensively. It enriched the study on Mars landing and the results provide theoretical support for the Chinese Mars exploration project. Moreover, the effectiveness of attitude control by thrusters is verified.

The rest of the paper is organized as follows.Section 2 presents the lander configuration and the arrangement of attitude control thrusters. In Section 3, the dynamics model is derived considering plastic deformation parts and nonlinear contact forces. Section 4 introduces the terrain generation method,and the dynamics model is verified by experiments in Section 5.In Section 6, Monte Carlo simulation in uncertain terrain is carried out and the effectiveness of the attitude control strategy is further analyzed. Finally, conclusions are drawn in Section 7.

2. Lander configuration

The Mars lander discussed in this paper is equipped with four sets of inverted-triangle landing gears, which are centrosymmetric and distributed around the main body. Each landing gear is mainly composed of a buffer rod, a primary strut,two secondary struts, and a footpad, as shown in Fig. 1. The top end of the buffer rod is fixed on the main body, and the other end is connected with the primary strut through a universal joint. One part of the secondary strut (outer cylinder) is connected with the main body by a universal joint, and the other part (inner cylinder) is connected with the footpad through a ball hinge. Due to the large cushioning capacity of the secondary strut, the primary strut does not provide cushioning and is fixed with the footpad.

As shown in Fig.2,a main engine is installed at the bottom of the lander for braking. When the footpad first contacts the Martian surface, the main engine begins to shut down. This landing strategy is called shutdown-at-touchdown mode. A certain amount of time is required for signal acquisition, decision, and shutdown operations after being touched down,which delays the shutdown process. In this case, the thrust of the main engine adversely affects the landing stability27.An attitude control system is applied to the lander to improve the landing adaptability and stability.10 thrusters are installed to control the attitude during landing, two of which are fixed on the upper plane of the main body, labeled as A1and A2,respectively. And the rest of the thrusters are equipped on the bottom of the main body, labeled from B1to B8. O-XYZ is the ground coordinate system whose OX axis points upwards against the gravitational direction, and the OZ axis points northwards of Mars.Om-XmYmZmis a body coordinate system fixed on the main body, whose OmXmaxis is the symmetry axis of the lander, and the OmYmaxis is parallel to the nozzle of A1. Oc-XcYcZcis a control coordinate system whose Ocand OcXcaxis coincide with the body coordinate system. The angle between the OcZcand the OmZmis 45 °.

Fig. 1 Landing gear configuration.

Fig. 2 Coordinate systems and thrusters’ position.

3. Dynamics modeling

3.1. Landing conditions

From a macro point of view, landing is a collision process between two objects. The landing conditions include relative velocity, relative angular velocity, and relative attitude. In order not to lose generality,landing conditions are divides into three groups.

(1) Slope parameters {αs, βs}.

Considering that Mars has a gravity of about 3.72 m/s2,the gravity effect is more evident than landing on the Moon(about 1.63 m/s2)or a small celestial body(almost no gravity).Therefore, the overall slope of the landing surface has a significant effect on the landing stability, and especially, it increases the risk of overturning. Here, the slope of the datum plane of the landing terrain αsand the orientation of the slope in OXYZ βs(as shown in Fig. 3) are considered. Since the surface of any celestial body is not entirely flat, these two parameters are only the overall slope of the landing area.In the dynamics model, rocks and craters are included.

Fig. 3 Landing conditions.

(2) Velocities {vh, βv, vx}.

Before contact,all parts of the lander have the same velocity relative to O-XYZ. Here, the horizontal velocity vh, the direction βv, and the vertical velocity vxare considered.

(3) Lander attitude{θZ,θX,θY}and angular velocities{ωX,ωY, ωZ}.

Rotate O-XYZ to make it parallel to Oc-XcYcZcby Z-X-Y rotation sequence and denote the rotation angles as θZ,θX,θY,respectively. The angular velocity vectors projected on Xc, Yc,Zcare denoted as ωX, ωY, ωZ, respectively. These attitude parameters are also the input of subsequent attitude control.

3.2. Nonlinear buffering force

The cushion material that connects the inner and outer cylinder of the secondary struts is the main energy-absorbing component through unidirectional plastic deformation. Its mechanical properties have been reported in detail in the literature17.Since we focus on the macroscopic dynamic characteristics of the lander,the detailed parameters and characteristics of this material are not repeatedly discussed here. The nonlinear relationship between the axial buffering force and the buffering stroke obtained in practice is shown in Fig. 4. Discretely measured points are interpolated by Cubic Spline Interpolation to simulate the force in the multi-body dynamics model. It should be noted that the buffering force is only capable when compression and the expression of the buffering force is

Fig. 4 Secondary strut buffering force.

where FSis the axial buffering force,DSis the buffering stroke,fiis the piecewise cubic polynomial.

3.3. Pseudo-rigid-body model of buffer rod

The buffer rod is used to absorb the impact between the primary strut and the main body. Plastic deformation occurs when the buffer rod is subjected to impact,and its deformation form is mainly bending. In order to simplify the dynamics model, the buffer rod is equivalently substituted with a pseudo-rigid-body model25,28–29,as shown in Fig.5.The buffer rod is modeled with two rigid links,a revolute joint,and a prismatic joint. Angle αEis used as the independent variable. The changes of moment MEand axial displacement SEwith respect to angle αEare obtained by experiment, as shown in Fig. 6.Analogously, discretely measured points are interpolated by Cubic Spline Interpolation.

3.4. Contact force

The contact force between the footpad and the landing terrain contains two parts, namely, the normal force and the tangential force. We adopt the nonlinear spring-damping model20to simulate the normal force, which can simulate the fact that the stiffness increases when the penetration depth increases.For the tangential force, some studies used the Coulomb friction for simulation.However,in the actual landing process,as shown in Fig.7,the footpad may get trapped in the soil.In this situation, in addition to the friction proportional to the positive pressure,there will be additional damping relate to tangential velocity. The coefficient of this damping will also increase with the penetration depth.

Since the traditional Coulomb friction model is no longer effective in accurately simulating the tangential mechanical characteristics, we modified the classic expression by adding a new depth-dependent damping term in the tangential force model as

Fig. 5 Pseudo-rigid-body model.

Fig. 6 Mechanical behavior of buffer rob.

Fig. 7 Contact model.

where Fnis the normal force that perpendicular to the landing terrain, K is stiffness coefficient, δ is the normal penetration depth, ε is exponent parameter, Ftis the tangential force that is opposite to the tangential velocity vt,μ is friction coefficient,Cnand ctare damping coefficients.The values of these parameters are shown in Table 1.

3.5. Thrust modeling

The attitude control thrusters are used to reduce the risk of overturning, which work for two seconds when the main engine has shut down. The control strategy concerns the instant attitude and angular velocity of the lander. When theydeviate from the ideal situation(θi=ωi=0,i=X,Y,Z),the control torques are computed as

Table 1 Parameters of contact model.

where a and b are control coefficients,kiis the angular velocity damping coefficient.

In practice, there are two situations to be considered. That is, there is a maximum thrust limitation for each thruster, and thrusters can only provide unidirectional force. Therefore, the thrust of a specific thruster is calculated according to the torque symbol and force arm, as

where FA1,FA2,FBiare the thrusts of A1,A2,and Bithrusters,respectively, FM1and FM2are the maximum allowable thrusts of A and B series control thrusters, and L1and L2are force arms. All control parameters are presented in Table 2.

3.6. Dynamics model

In this paper, the Lagrange dynamic equations are used to deduce the dynamics model of the Mars lander. An overview of the three-dimensional lander simulation model is shown in Fig. 8. After simplification, the whole model contains 29 rigid parts and 40 groups of constraints, including 12 translational joints, 12 universal joints, 8 spherical joints, 4 revolute joints,and 4 αE-SEconstraints. Some necessary parameters of the lander are presented in Table 3.

The centroid position and Euler angle of each rigid part with respect to O-XYZ are considered generalized coordinates,and the dynamics equations are expressed as

Table 2 Parameters of control coefficients.

Fig. 8 Overview of three-dimensional lander simulation model.

Table 3 Parameters of the Mars lander.

where M is the mass matrix, C is 156 × 1 constraint function vector, q is 174 × 1 generalized coordinates vector, Cqis the first partial differentiation of C with respect to q,λ is the vector of Lagrange multipliers,Qvis the vector of generalized inertial forces, the subscript t denotes partial differentiation with respect to time, and Qeis generalized external force vector, as

where G is gravity,FSis buffering forces of the secondary strut,MEis the equivalent moments of buffer rods, FCis contact forces, and FTis thrust.

The Runge-Kutta method is used to solve the dynamic differential equations.The integral step is 0.001 s with a total simulation time of 10 seconds.

4. Uncertain landing terrain modeling

In order to comprehensively evaluate the landing performance of the lander landing in irregular terrain, we aim to obtain a local terrain model which can be used in dynamics simulation.It is rather difficult to express the terrain surface with an analytical formula.However,according to the existing techniques,the terrain can be modeled with two common characteristics:meteor craters and rocks.Since the comprehensive terrain data of Mars is not available,the terrain model is established,referring to the statistical data of the distribution of lunar rocks and craters. According to the literature30, typical distributions of craters and rocks on the lunar surface can be approximately expressed as

where Dcis the diameter of a crater, Dbis the diameter of a rock, Ncis the accumulated number of craters larger than Dcper 1 m2, Nbis the accumulated number of rocks larger than Dbper 1 m2. Craters with Dc< 0.2 m and rocks with Db<0.05 m are removed because tiny craters and rocks have little effect on the process of landing.

The shape of a crater is usually irregular,which brings great difficulties in modeling. In this paper, we assume that craters are all in the regular shape expressed as

Examples of a crater and rock are shown in Fig. 9.

Fig. 9 Examples of crater and rock.

Fig. 10 Generated landing terrain.

Now that we have the expressions in Eqs. (8)-(12), craters and rocks are randomly generated on a 50 m × 50 m landing plane using the proposed Algorithm 1 in the following. And the generated terrain is shown in Fig. 10.

Algorithm 1 Randomly terrain generation 1: Inputs: Given terrain length l1 and width l2 2: Step1: Compute the set of craters diameters Ιc 3: Ιc =?4: Terrain area S = l1 l2 5: Solve equation:Nc(Dc)S=1,obtain Dcmax 6: while Dcmax > 0.2 7: Generate random number p ?[0,1]8: Dnow = pDcmax 9: Ιused = {Dc | Dc ?Ιc, Dc > Dnow}10: Nused = card (Ιused)// card(·) return the number of elements in a set 11: Nt = Nc(Dnow)12: Nr = Nt- Nused 13: if Nr < 0 14: Dcmax = Dnow 15: Dnow = 0 16: else 17: Ιc ←Ιc∪{Dnow}18: end 19: end 20: Step2: Compute the set of rocks diameters Ιb 21: ... (Similar to step1)22: Step3: Generate random position x0 ?[0, l 1] and y0 ?[0, l2] for all elements of Ιc and Ιb 23: Step4: Calculate the shape for all elements by Eqs. (10)-(12)24: Step5: Superimpose and output terrain

5. Model validation

Model updating technology31–33is used to improve the accuracy of the dynamics model. In order to reduce the difficulty of the experiment,the landing terrain is a plane,and the effects of the thrusters are not included in the experiment.Experiment facilities are illustrated in Fig. 11.

The experiment is conducted with the following steps:

(a) Connect the prototype with the double-rocker linkage.Pull it to a certain height and then release the lander.

(b) Unchain the prototype from the double-rocker linkage at the lowest position to give the prototype an initial horizontal velocity.

(c) Prototype lands on the slope.

Fig. 11 Experiment facilities.

Landing conditions are simulated by adjusting the attitude control device, the release height of the prototype, and the slope of the landing terrain. A high-speed camera system is used to measure the landing velocity and attitude.Meanwhile,the minimum distance between Omand the landing terrain(HMHM) and the secondary struts buffering strokes (DSDS)are measured. More details are discussed in the literature34.Two sets of experimental data were used for validation, and landing conditions are presented in Table 4.

In case 1,there is no slope(0°)and the landing speed is relatively low, so the landing performance can be preliminarily tested. Case 2 increases the slope (8°) and the landing speed,which can further test the performance under adverse conditions.

Landing parameters in the simulation are set according to Table 4. The comparison between simulation and experiment is shown in Fig. 12. The simulation results show that the proportion trend of DSbetween different legs is close to that of the experiment.The average error of DSis 2.80 mm,and the average relative error is 10.5%. The average error of HMis 2.17 mm with a relative error of 0.25%. For the experiments,the measurement accuracy and manufacturing error can affect the results. Therefore, since relative errors are within the acceptable ranges, the accuracy of the dynamics model is verified.

6. Stability analysis

6.1. Stable criterion

This paper focuses on the stability of soft landing under the shutdown-at-touchdown mode, which mainly considers antioverturning ability, buffering performance, and the attitude after landing.Observations indicating the landing performance and stability are chosen as follows.

(1) Lmin(>50 mm). Minimum distance from the center of mass to the overturning plane12. The overturning plane is defined as the vertical plane passing through the centers of two adjacent footpads.This variable indicates the ability to resist overturning. When Lminis close to zero,the lander has a high risk of overturning.

Fig. 12 Simulation and experiment results.

(2) Hmin(>200 mm). Clearance between the nozzle of the main engine and the landing terrain. If Hminis small,the nozzle may collide with rocks, resulting in damage.

(3) αEmax(<30°).The maximum rotation angle of buffering rod.

(4) DSmax(<100 mm). The maximum buffering stroke of secondary struts.

(5) Amax(<150 m/s2). The maximum acceleration of the main body.

(6) ψ (<20°). The angle between the Ym-Zmplane and the slope after landing.If ψ is large, when the rover departs from the lander,the rover may collapse during the transfer process and result in mission failure.

6.2. Monte Carlo simulation

Due to uncertain factors such as external disturbance, landing velocities and attitude cannot be accurately determined. Generally, these parameters are constrained in particular ranges.Therefore,the simulation under certain conditions is not effective to completely evaluate the landing stability, especially when effects of random terrain are considered. Monte Carlo simulation35is an alternative to deal with such uncertain problems. This method is based on statistical sampling for numerical calculation,and the more simulation samples are,the more accurate the results are.The computational cost of multi-body dynamics simulation is relatively low and suitable for a large amount of computations.

Fig. 13 The Monte Carlo Simulation process.

The Monte Carlo Simulation analysis process is shown in Fig. 13. The first step is to generate the 50 m × 50 m landing terrain using Algorithm 1.Then the next step is to generate the landing conditions and select the landing position randomly on the landing terrain with equal probability. According to the actual design capability,the probability distribution of landing attitude and landing velocity are listed in Table 5. Here, βvobeys uniform distributions, and other landing conditions are assumed to obey normal distributions.

It is noteworthy that when the lander lands, each footpad only contacts a small area of the landing terrain. Therefore, it is necessary to obtain local slope,craters,and rocks information according to the landing site before dynamic simulation. The final step is to carry out dynamics simulation and output results.In this paper,a total of 1500 simulation samples are obtained.

The distribution histograms of the Monte Carlo Simulation results are shown in Fig.14.As shown in Fig.14(a),in 5.0%of the simulation samples,the lander overturns(Lmin<50 mm).Most Hminis distributed in the safe range of[300,500]mm,as shown in Fig. 14(b), and 0.6% of Hmin< 200 mm because of rocks.The maximum rotation angle of the buffering rod αEmax,the maximum buffering stroke of the secondary strut DSmaxand the maximum acceleration of main body Amaxare all in the safe range and far away from the dangerous threshold,as shown in Fig. 14(c)-(e). In some cases, the angle between the YmOmZmplane and the landing terrain ψ exceeds the safe threshold,with a ratio of 0.5%(overturning samples are excluded).

The simulation results suggest that the lander has a stable landing ratio of 94.5% when landing randomly on the generated landing terrain. The overall insecurity rate is 5.5%, and the primary failure forms are overturning and abnormal posture caused by rough topography, while the buffering design is redundant.

Fig.15(a)shows all the landing sites on the generated terrain and Fig.15(b)shows regional stability with a 5 m×5 m division.Overall,most areas have 100%stability and the unstable landing sites can be divided into two types. (1) Landing sites are located at the inner edges of large craters (craters A and C),where a high equivalent slope exists.The regional stability of this part is about 80%.(2)Most of the unstable landing sites are near the area between crater B and C,where the terrain is rugged,and the heights vary dramatically.The regional stability of this area decreased significantly to less than 50%.

6.3. Effects of attitude control

In previous analysis results, Lminis the most critical indicator of landing stability.The primary purpose of attitude control is to improve the anti-overturning ability. In order to quantitatively evaluate the effect of attitude control thrust on stability,another 1500 dynamics simulations are performed. The landing parameters of these simulations, including attitude, velocity, and landing sites, are the same as those in Section 6.2 that the only difference is that attitude control thrusters are deactivated.

Table 5 Probability distribution of landing conditions.

Fig. 14 Histogram of Monte Carlo Simulation outputs.

Fig. 15 Landing sites on generated terrain and regional stability.

The simulation results are compared with the previous. As expected, when deactivating the attitude control thrusters, the overall stability is decreased with the overturning ratio increasing from 5.0%to 8.9%,as shown in Table 6.The stable ratio is 94.5% when the attitude control thrusters are used, and it decreases to 90.7% when thrusters are deactivated.

Concretely,the box plots, which contain the median,quartiles, extreme values, and outliers, are shown in Fig. 16. Most outputs are worse except for Hminwhen the attitude control thrusters are closed. Take the median as an example: LminLmindecreased from 852.2 mm to 838.9 mm, αEmaxαEmaxincreased from 3.86° to 3.87°, DSmaxDSmaxincreased from 14.6 mm to 14.8 mm, AmaxAmaxincreased from 43.7 m/s2to 44.3 m/s2, and ψ increased from 1.50° to 1.63°. The variation trend of the quartile, the maximum value, and the minimum value are consistent with the median.

Overall,when the thrusters are used,the landing stability is improved by 3.8% in the whole 50 m × 50 m landing area.Referring to Eq.(3), when the landing surface is relatively flat,the thrust is usually low, and in this situation, the differencebetween landing with thrust on and off is slight. When the landing surface is rough, the moment provided by thrusters can significantly enhance the anti-overturning ability, even if the thrust action time is short.Fig.17 shows the unstable landing sites when thrusters are on and off. In a bad local area around the crater,it can be observed that the unstable samples are significantly reduced when using thrusters.The comparison results validate the positive effects of attitude control thrusters.

Table 6 Stability when thrusters are on and off.

Fig. 17 Unstable landing sites when thrusters are on and off.

Fig. 16 Box plot of simulation outputs with and without control thrusters.

7. Conclusions

(1) A multi-body dynamics modeling method for a Mars lander is developed in this paper. The pseudo-rigid body method is used to model the plastic deformation parts for fast computation. And the nonlinear springdamping model and modified Coulomb friction model with correction terms are used to model the contact. In addition,10 GNC attitude control thrusters are activated to improve the anti-overturning ability during landing.

(2) The effectiveness and accuracy of the dynamic modeling method are verified by experiments.The average error of secondary struts buffering strokes is 10%,and the average error of minimum distance between the origin and the landing terrain is 0.25%.

(3) A random landing terrain generation method based on statistical data is proposed.A 50 m×50 m random terrain is obtained by linear superposition of randomly generated craters and rocks with regular shape. The establishment of random terrain benefits the analysis of the terrain adaptability of the Mars lander.

(4) 1500 sets of soft landing simulations are performed to analyze the landing stability and the effects of attitude control thrusters under uncertain terrain. The results show that the insecurity rate is 5.5%, in which the primary failure forms are overturning and abnormal posture caused by rough terrain. Another 1500 simulations without thrusters are performed. And the results suggest that the attitude thrusters can improve the landing performance under rough terrain.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This study was supported by the Joint Fund of Advanced Aerospace Manufacturing Technology Research(No.U2037602).