Influence of the Lander Size and Shape on the Ballistic Landing Motion

Abstract

This article investigates the ballistic landing motion and final distribution of the landers in different sizes or shapes near the small celestial body. Three typical shapes, including cubic, cuboid, and cylindrical, are considered for the landers deployed to a tri‐axial ellipsoid model. The Polygonal Contact Model (PCM) is used to detect the contact/collision, where the Hertz model is applied to calculate the continuous contact force. Different‐sized cubic landers (in the edge length of 20, 30, 40, and 50 cm) are numerically simulated to examine how the lander size influences its dynamics. The landing motion of the cuboid‐ and the cylinder‐shaped landers are then analyzed in the same technique. The heights of these asymmetrical landers are assumed to be 25, 30, and 35 cm, respectively, to illustrate the shape effect. Monte Carlo simulations are implemented for various landers to account for the surface motion randomness. The final dispersion, the outgoing velocity after the collision, the horizontal transfer distance, and the settling time are taken to be critical indicators for discussing the landing behavior, which can provide implications for the probe design of future missions.

Introduction

Small celestial bodies are the exploration targets for many deep space missions (Biele et al., 2015; Kawaguchi et al., 2008; Li et al., 2013; Watanabe et al., 2017; Yeomans et al., 2000). Small body missions can help humans discover the formation mysteries of the solar system (Boehnhardt et al., 2017), exploit new storage for mineral products or organics (Nadoushan et al., 2020), and test the planetary defense technologies designed to protect the Earth’s environment (Cheng et al., 2016). Although the remote orbital observations have provided some beneficial external characteristics, the surface environment, and the inner structures of small bodies are unknown except for implementing the in-situ surface exploration. The hopping probes have been demonstrated to be feasible instruments for in-situ exploration (Ulamec et al., 2011). For example, the rovers, named MINERVA-II-1A, MINERVA-1B, and the lander Mobile Small body Surface Scout (MASCOT) have been successfully deployed by Hayabusa 2 (Scholten et al., 2019; Soldini et al., 2020; Yoshikawa et al., 2020). They landed on the surfaces of the small body 162173 Ryugu, sent back the pictures and videos of the rocky surfaces, examined the surface mechanical properties, and implemented a series of scientific missions, such as investigating the multispectral reflectance, the thermal characteristics, and the magnetic properties (Ho et al., 2021). The release altitude above the Ryugu surface for the two probes is estimated at 40–50 m, respectively (Yoshikawa et al., 2020), which is far less than the orbital observation distance. The reason for the lower deployment height is that the irregular micro-gravity near the small body makes the probe very easy to escape without control. The rover MINERVA of Hayabusa is an example that escaped from the small body 25143 Itokawa because the spacecraft deployed the probe prematurely under false instruction (Yoshimitsu et al., 2006). Additionally, the small body surfaces are rough and complicated, with many irregular rocks and weathering debris, which is also an obstruction to the success of the surface mission (Lauretta et al., 2019; Mahaney & Kapran, 2009; Michikami et al., 2019).

In this way, parametric simulations of the landing motion before carrying out a realistic mission are essential. In some earlier studies, the lander was usually regarded as a test particle for discussing its trajectory (Bellerose & Scheeres, 2008; Çelik & Sánchez, 2017; Ferrari & Lavagna, 2018; Song et al., 2020; Y. Zhang et al., 2019). Hockman and Pavone (2020) contributed to the landers’ motion planning based on the Lambert problem. Yang and Baoyin (2015) contributed to the fuel-optimal control problem for the soft landing on an irregular small body and established an improved homotopic method. Wen et al. (2020) proposed the hop reachable domain of the particle on Bennu and comet 67P/Churyumov-Gerasimenko. To provide high-fidelity simulations of the rigid-body lander, Tardivel et al. (2014) first modeled the spherical probe hopping on the polyhedron model of Itokawa. As a follow-up, Van Wal et al. (2017, 2020) continuously examined the influence of the lander mass distribution and the probe shape on the settling time and the surface dispersion. Five platonic solids were discussed and compared to find how the lander shape affects the deployment dynamics (Van Wal et al., 2020). However, all shapes are given the same density and mass, and only one size is selected and discussed for each lander shape. The interaction between the nonspherical lander and the regolith has also been studied by applying the soft-sphere model (Çelik et al., 2019; Maurel et al., 2018; Thuillet et al., 2018). In addition to these parametric simulations, the mechanics of the spacecraft-regolith interactions have been tested in experiments (Murdoch et al., 2021; Van Wal et al., 2021). Van Wal et al. (2021) tested Stronge’s impact model using a rectangular assembly with air bearing footpads that was placed on an inclined granite table, to create an artificial, two-dimensional microgravity environment. Murdoch et al. (2021) performed the low-velocity collisions into different types of granular materials in the terrestrial gravity and the reduced gravity. The experimental data support a drag model for the impact dynamics composed of both the hydrodynamic drag force and the quasi-static resistance force, indicating that regolith has a more fluid-like behavior in low gravity.

This article aims to analyze the influence of the rigid lander size and shape on the landing motion. The lander is taken to be a rigid body where the surfaces of small bodies are treated as elastic layers by assuming the contact areas are in rocky terrains. A homogeneous ellipsoid is adopted as a general model approximating irregular small bodies, whose parameters are taken resembling the comet 133P/Elst-Pizarro. This comet has significant scientific values, which was chosen to be the target of the Chinese asteroid exploration mission in 2019 (X. J. Zhang et al., 2019).The properties and exploration significance of this comet will be introduced in Section 3.1 in more detail. Three different shapes of landers, cube, cuboid, and cylinder, as shown in Figure 1, will be discussed. As for the cubic landers, with the same mass and the moment of inertia, their edge lengths are changed from 20, 30, 40, to 50 cm, respectively, for analyzing the individual influence of the lander size. Afterward, keeping the mass identical, the 20 cm-sized cubic lander is transformed into the cuboid lander to study the shape effect. One edge of the 20 cm-sized cube is enlarged to 25, 30, and 35 cm, respectively. The cuboid landers are assumed to be homogeneous when calculating the moment of inertia. Keeping the corresponding mass and the moment of inertia identical, the lander shape is transferred from the cuboid to the cylinder with uniformly distributed mass by resolving the radius and height inversely. The Monte Carlo simulations with 1,000 trajectories are made for these 10 different landers due to the stochastic motion near the small body. The final dispersion, settling time, horizontal transfer distance, and outgoing velocity direction are the critical parameters when analyzing the influence mechanism.

Figure 1. Ten typical landers in different sizes and shapes.

This article is organized as follows. Section 2 introduces the fundamental dynamics used to describe the coupled orbit-attitude landing motion. The PCM method for determining the collision and the Hertz model for calculating the continuous contact force are also briefly introduced for the integrity of the paper. Section 3 highlights the parametric simulations of the lander with different sizes and shapes, where system parameters of the ellipsoid are chosen according to the comet 133P. The specific influence of the lander shape and size on the trajectory evolution is summarized and compared. The coupled influence of the lander mass and size is analyzed. Additionally, the landing mechanisms are discussed according to the numerical results. Finally, the conclusion is summarized in Section 4.

Dynamical Equation Coupling Orbit and Attitude Equations

In this study, the coupled orbit and attitude model of the rigid lander is analyzed, that is, the velocity and rotation rate of the lander are changed simultaneously. According to the absence or presence of contact with a small body, the landing motion can be classified into two modes: the ballistic transfer mode and the collision mode. The gravitational force is supposed to act on the centroid of the lander. In this case, the gravity-gradient torque is zero, and the contact torque significantly alters the lander’s rotation rate during the collision. To discuss the lander’s attitude, the small body mass center principle frame O-XYZ and the lander rotating frame p-xyz are defined in Figure 2, respectively. The sizes of the small celestial body and the landers are not scaled. The origin of the small body principle axis frame O-XYZ is at the centroid of the small body, with three orthogonal axes aligning the principal axes of inertia. The parameter ω_α represents the angular velocity of the small body relative to the J2,000 inertial frame. In the short period of the landing motion near the small body, the angular velocity ω_α is regarded as a constant parameter. The lander rotating frame p-xyz has its origin at the centroid. The three axes are along with the maximum moment of inertia, the intermediate moment of inertia, and the smallest moment of inertia, respectively. The parameter v is the velocity of the lander’s center of mass and ω_β is the angular velocity of the lander.

Figure 2. Coordinates of the small body and the lander.

When studying the landing motion, the micro-gravity of the small body, the centrifugal force, the contact force, and the tangential friction force are considered in this paper. The solar radiation pressure force and the solar or planet tides are omitted in the equations due to the relatively short timescale of landings. In the ballistic transfer, 4 vectors with 13-dimensional variables represent the lander state in Equation 1 [Image Omitted. See PDF]

Here, the variables r and v are the position and the velocity of the lander’s center of mass relative to the small body frame, respectively. The angular velocity ω_β of the lander and the quaternion q are applied to describe the lander’s attitude motion. Based on position r, the dynamical equation of the lander is expressed in Equation 2 [Image Omitted. See PDF]

Here, the parameter $$-\nabla U$$is the gravitational acceleration derived from the polyhedron method. In order to build a precise model suitable for different kinds of asteroids with complex shapes, the polygonal model is utilized for calculating the weak and irregular gravity model of the asteroid. The parameters G and ρ represent the gravitational constant and the density of the small body, respectively. The variables r_e and r_f are the vectors from the field point to the polyhedral points on the edge e and the face f, respectively. The parameters E_e and F_f are the constant edge dyad and face dyad, and L_e and θ_f are the edge factor and the trigonal solid angle. The subscript “Edge” and “Face” are the sets of all edges and facets of the polyhedron model (Werner & Scheeres, 1997). The components of the above vectors are expressed in the small body frame.

Furthermore, the units of the variables are uniformed in Equation 3 for the convenience of comparing different landers subsequently. The time unit [T] and the distance unit [L] are the rotation period and the equivalent radius of the small body. [Image Omitted. See PDF]where the superscripts “∼” and “∧” represent the normalized vector and the unit vector, respectively. The parameters T_a and V_a are the rotation period time and the volume of the small body. Additionally, the introduced parameter κ describes the integrated effect of the gravity field and the centrifugal force of the small body. When analyzing the lander attitude, the angular velocity ω_β and the quaternion q are applied, where the attitude dynamical equations are expressed in the lander frame p-xyz in Equation 4 [Image Omitted. See PDF]

Here, [I] is the moment of inertia of the lander. The parameter L_g is the torque of the gravitational force, which equals zero in this study, with the gravity going through the centroid of the lander.

In the collision mode, the lander’s motion is mainly affected by the contact force that comprises the normal support force and the tangential friction. With the assumption of the elastic layers on the surface of the small body, the Hertz model is adopted to calculate the continuous contact force during the collision (Lankarani & Nikravesh, 1994). The elastic layers guarantee that the deformation of the asteroid surfaces can recover after the contact. The total collision period can be divided into the compression and the restitution phases. According to this model, two parameters K and D are introduced before calculating the normal force in Equation 5 [Image Omitted. See PDF]where K is the stiffness of the surface, the parameters γ_i, Y_i, and R_i (i = 1, 2) are the Poisson ratio, the Young’s modulus, and the radius of curvature of the two contact bodies, respectively. Usually, the value of the surface stiffness K can be directly tested in the experiment or derived from the actual rebounding trajectory (Yano et al., 2006). As a preliminary study, the surface stiffness is selected as a constant value (70 N·m^−1.5). The parameter δ is the penetration depth, $${\dot{\boldsymbol{\delta }}}_{(-)}$$ is the relative normal velocity before contact, and the restitution coefficient e is a parameter between 0 and 1. The normal contact force is yielded in Equation 6 [Image Omitted. See PDF]

Here, N represents the normal force, including two parts: $$K{\delta }^{\frac{\mathrm{3}}{\mathrm{2}}}$$ represents the elastic contact force, which is directly related to the amount of local indentation, $$D{\Vert}\dot{\boldsymbol{\delta }}{\Vert}$$ is applied to describe the energy dissipation in the form of internal damping of colliding solids. $$\hat{\boldsymbol{n}}$$ is a unit vector representing the normal direction, and $$\dot{\boldsymbol{\delta }}$$ is the normal penetration velocity. The critical parameter e is defined in Equation 7 [Image Omitted. See PDF]

In the above equation, $${\dot{\boldsymbol{\delta }}}_{(+)}$$ is the relative normal velocity after contact. The collision integration terminal can be controlled by e. When the normal magnitude of v exceeds $$e{\dot{\delta }}_{(-)}$$ during the restitution integration, the collision period reaches the terminal, and then the lander re-enters the ballistic transfer mode. Without the exact surface property, e is set as 0.5 temporarily (Jourdain, 1917).

The tangential force F_f is solved based on the relative tangential velocity v_t and the normal contact force in Equation 8 [Image Omitted. See PDF]

This equation is a regularized version of the Coulomb friction law. In Equation 8, μ is the friction coefficient, N is the magnitude of the normal force N in Equation 6, and V_d is the chosen velocity threshold. When the magnitude of v_t falls below V_d, the friction force is faded out quadratically to avoid the set-value state of the static friction (Hippmann, 2004). Finally, the contact force F_d is the summarization of the normal force and the tangential force: F_d = N + F_f. With the contact force, the relative derivative of X in the collision process is given in Equation 9 [Image Omitted. See PDF]where m represents the lander mass and L_d is the torque of the contact force F_d.

Detection of Contact and Settling

The Polygonal Contact Model (PCM) method is applied to detect the active collision areas of the polyhedral models (Hippmann, 2004). The bounding volume (BV) hierarchies are utilized to detect the collision, which can reduce the calculation burden to a great extent. The BV hierarchies are illustrated utilizing an ellipsoid in Figure 3. The BV hierarchies belong to a Binary tree data structure. The cuboids are applied as BVs, and the root cuboid includes the whole triangular facets of one polyhedral model. The root cuboid is recursively divided into two child BVs, each including half of the triangular facets. This operation is repeated until the leaf BV node only includes one triangular facet. With the BV hierarchies, the collision detection starts with the two root BV cuboids. If the root cuboids interact, the collision tests between their child BVs are carried out. Once the upper BV cuboid is excluded from the collision, it is not necessary to detect its child nodes. If the triangles in the leaf nodes have interacted, the points of interaction on the edges can be obtained. The intersection polygons are constructed by distance calculation of interaction points. Furthermore, the inner triangular facets surrounded by the intersection polygons are searched, which are the active collision areas for calculating the contact force.

Figure 3. Bounding volume hierarchies and the active contact facets.

When the lander leaves the surfaces after interacting with the small body, the next potential hopping height h_e is calculated in Equation 10 according to the outgoing velocity and the effective potential. When h_e is small enough, the lander is considered settled. [Image Omitted. See PDF]where V represents the effective potential of the small body, including the influence of gravitational and the centrifugal terms. The effective potential is a conserved quantity for the conserved dynamical motion of the lander in the small body frame, which can be added to the lander’s kinetic energy relative to the small body frame (Scheeres et al., 2016). A height threshold h_ε is given in the following simulations, which equals 1 cm. During each collision, the estimated hopping height h_e is calculated. If h_e > h_ε, the lander propagates continuously. Otherwise, the lander is assumed to settle down after this collision.

Parametric Simulation Model of the Central Body

The comet 133P/Elst-Pizarro is regarded as the simulation background in this article. It is the first discovered main-belt comets (MBCs), which has been once selected as the possible exploration target in several missions (Gao et al., 2021; Jones et al., 2018; Snodgrass et al., 2018). Exploring MBCs is of great significance for scientific exploration and may help explain the mystery of the origin of the terrestrial water on Earth (Hsieh and Jewitt, 2006). However, the accurate shape model of comet 133P is not available without close-range detection. Therefore, a homogeneous ellipsoid is applied to approximate comet 133P based on the obtained physical characteristics from the ground-based observation. The density of this comet is about 1,300 kg·m⁻³, and its equivalent radius equals 1.9 ± 0.3 km with three semi-axes are 2.3, 1.6, and 1.86 km, respectively (Hsieh et al., 2009). The polyhedral ellipsoid model applied in this article has 42,914 vertices and 85,824 triangular facets, with the mean edge length in 36 m. The ellipsoid model's effective potential V in Equation 10 is also calculated and illustrated based on the observed data in Figure 4. The infinite position from the central body has zero gravitational potential. In Figure 4, the effective potential is higher in the polar regions of the ellipsoid and lower in the equatorial region. If a lander is released near the polar region, it tends to move toward the equatorial region spontaneously (Y. Zhang et al., 2019).

Figure 4. Effective potential distribution of the central ellipsoid model.

Cubic Landers in Different Sizes

For discussing the influence of the different-sized cubic landers on the trajectory evolution, their edge lengths are chosen to be 20, 30, 40, and 50 cm, respectively. Without loss of generality, their mass and the moment of inertia are selected based on the lander MASCOT (Ho et al., 2017). Particularly, they are identical for the landers with different sizes. The lander’s mass equals 10 kg, and the moment of inertia is [0.067, 0.067, 0.067]^T kg·m². In this case, the 20 cm-sized lander has a homogeneous mass distribution in the highest mean density, while the larger-sized lander’s mass distribution becomes more concentrated at its geometric center. From a physical point of view, the change of the mass distribution can be achieved by altering the inner positions of the onboard instruments, which is out of scope of this study. The identical mass and inertia of the landers ensure the same initial kinetic energy and effective potential energy when the release situation is determined. The effective potential energy is the product of the lander mass and the effective potential V. Only the lander size is varied for analyzing its influence on the lander behavior.

The initial states (including the release position, velocity, quaternion, and angular velocity) of the simulations are listed in Table 1. Monte Carlo simulations are utilized in this study. In general, a great deal of randomness of the landing motion exists near the small body surfaces. It is hard to predict the dynamical propagation of the lander through only one specified trajectory (Van Wal et al., 2017). Therefore, the Monte Carlo simulation is necessary to increase the reliability of results when the behavior of different landers is discussed. According to the Monte Carlo method (Metropolis & Ulam, 1949), one thousand simulations are performed for each lander size with tiny perturbations of the lander’s initial position and velocity.

Table 1 Parameters in the Simulations

Parameters	Symbols	Values
Cubic lander edge	rc	0.20–0.50 m
Cubic lander mass	m	10 kg
Cubic lander moment of inertia	I	[0.067, 0.067, 0.067]T kg·m2
Coefficient of Hertz model	e	0.5
Coefficient of friction	μ	0.45
Initial position	x	[9.59, 1,883.81, 134.76]T m
Initial velocity	v	[0.00,0.00, 0.00]T m/s
Initial angular velocity	ωβ	[0.00,0.00, 0.00]T rad/s
Initial quaternion	q	[1, 0, 0, 0]T
Initial position uncertainly (standard deviation)	xσ	[1, 1, 1]T m
Initial velocity uncertainly (standard deviation)	vσ	[0.0002, 0.0002, 0.0002]T m/s

The cubic landers are released from the identical initial position with zero velocity and no spin. The positions and velocities refer to the state of the centroid hereafter. The release altitude in this study is about 20 m above the equatorial surfaces of the asteroid, which is lower than the previous studies. In previous studies, the landers were numerically released from hundreds of meters above the surfaces (Tardivel et al., 2014). In these cases, the significant change of the effective potential plays a coupling role in the landing motion propagation, resulting in very diffuse trajectory settling positions. For example, the settling positions of landers released from [0, 0, 200]^T m can cover more than 13% of Itokawa’s total surfaces (Van Wal et al., 2020). To reduce the influence of the effective potential change, the lower release altitude is utilized in this study to focus mainly on the lander shape’s effect. Additionally, referring to the previous articles, the initial velocity of the lander is set as zero (Van Wal et al., 2017; Zhang et al., 2021).

Taking the 20 cm-sized and the 30 cm-sized cubic landers as an example, the trajectory evolutions of the 1,000 trajectories are illustrated in Figure 5 with some representative points and trajectories, respectively. The settling positions are located in a local hexagonal region constituted by six triangular facets. The total area of this region is about 0.013 km², and the latitude and the longitude are [89.1°, 91.3°] and [2.8°, 4.7°], respectively. According to Section 3.1, the superficial area of the central body is around 46 km². Therefore, the effective potential influence on the settling position can be reasonably neglected. Lander size is regarded as the only factor determining the energy loss, the trajectory evolution, and the final dispersion.

Figure 5. Trajectory evolutions of the cubic landers. (a) A 20 cm-sized lander. (b) A 30 cm-sized lander.

In Figure 5, the green asterisks, the blue dots, the yellow crosses, and the red curves (color online) represent the initial points, the settling positions, the first landing points (referred to as TD1 points hereafter), and the hopping trajectories, respectively. The positions of TD1 points are identical for different landers, and these points are in the region with the longitude and latitude scopes are [89.69°, 89.88°] and [3.85°, 4.06°], respectively. The lander with a smaller edge length tends to hop higher during the transfer process. The 20 cm-sized cubic lander has the most variable settling positions, and some of the trajectories hop upward from the low latitude to the high latitude. On the contrary, when the lander edge increases from 30 to 50 cm, all the trajectories transfer toward the lower latitudes. This difference reflects that the attitude and velocity of the lander are altered significantly by the size change during the collision, and the outgoing result after contact is sensitive to the state before the collision.

The coupled orbit-attitude of one representative trajectory of the 30 cm-sized cubic lander is illustrated in Figure 6. The white trajectory is enlarged in Figure 6a, focusing on the rebounding motion before the settling. The change of the translational energy and the rotational energy of the lander are illustrated in Figure 6b, respectively. The initial translational and rotational energies are zero because the lander is released with zero velocity and no spin. During the free landing before the first collision, the translational energy increases gradually to 0.043 J. When the lander collides with the surfaces, the velocity decreases in the compression period and increases partly in the restitution stage. Therefore, the translational energy decreases after each collision and eventually becomes zero. The rotational rate remains zero before the first collision because the gravity-gradient torque is zero. The conversion between the translational energy and the rotational energy is observed at the first collision, where the rotational energy increases from zero to 0.0014 J. Subsequently, the spin rate of the lander decreases due to the contact force. It turns into zero at the settling position.

Figure 6. Trajectory evolution of the 30 cm-sized cubic lander. (a) Trajectory evolution of the lander. (b) Translational energy and rotational energy of the lander.

The settling positions’ longitudes and latitudes are calculated and illustrated to analyze the concentration degree in Figure 7. The dark red dots (color online) represent the settling positions of the 20 cm-sized landers. The settling positions of the other three landers are omitted for brevity. Taking the longitude and latitude as two variables, the 95% confidence ellipse describes the concentration degree. With a given confidence coefficient, the error ellipse is determined by the mean values and the covariance matrix of the settling points’ longitudes and latitudes (Zhu et al., 2019). In the Monte Carlo simulation, the 95% error ellipse indicates that the probability that one trajectory settling position will fall in this ellipse is 95%. The geometrical area of the 95% error ellipse is calculated based on the coordinate values of the longitudes and latitudes directly. Although the ellipse area does not have a specific physical meaning, it can reflect the 95% confidence interval distribution. The smaller the 95% error ellipse area, the more concentrated the settling positions, reflecting that this lander is less sensitive to the initial release disturbance. Settling positions of the 20 cm-sized landers in Figure 7 distribute in the region with the longitudes and latitudes between [89.70°, 90.10°] and [3.80°, 4.20°], respectively. When the lander edge increases from 30 to 50 cm, the settling positions have similar longitudes that change from 89.90° to 90.20°, while the latitude scopes are approximately [3.55°, 3.77°], [3.60°, 3.80°], and [3.62°, 3.82°], respectively. Latitudes of the final distribution areas gradually move away from the equator. Additionally, the 95% confidence ellipses' areas of the landers with 20–50 cm edge equal to 0.066, 0.027, 0.025, and 0.032, respectively. The area of the error ellipse can reveal the concentration degree of the final distribution, and the 20 cm lander has the most divergent settling positions. The 30, 40, and 50 cm landers have similar concentration degrees, and the trajectories of the 40 cm lander are slightly more gathered.

Figure 7. Final distributions of different-sized cubic landers released from 20 m.

In order to explain the differences in the settling position distributions, the outgoing velocities at TD1 points, the settling times, and the transfer distances are illustrated and compared. Because of the lower release altitude, the outgoing velocities of the different-sized cubic landers have similar magnitudes. However, the directions of these velocities are significantly different. The directions of the outgoing velocities after the first collisions process are described based on the two parameters ϕ and λ in Figure 8, which are referred to as the elevation angle and the azimuth angle.

Figure 8. Two angles for describing the outgoing velocity direction.

The local coordinate Q−ξ₁ξ₂ξ₃ is constructed by two edge vectors and the normal vector of the triangular facet where the landers landed. The three facet vertices of the first landing triangle are supposed to be a, b, and c in the counterclockwise order. The origin Q of the local coordinate is set as the centroid of the triangular facet. First, the normal vector ξ₃ of this facet is calculated in Equation 11 [Image Omitted. See PDF]

Based on ξ₃, ξ₁, and ξ₂ are defined as Equation 12 [Image Omitted. See PDF]

In the local facet coordinate, a given velocity v can be solved in Equation 13 [Image Omitted. See PDF]

Therefore, given the velocity v, and ξ₁, ξ₂, and ξ₃, the angles ϕ and λ can be computed inversely from Equation 13. The angles ϕ and λ of the outgoing velocities at TD1 points are illustrated for different-sized cubic landers in Figure 9. The Kernel Density Estimation (KDE) toolkit in Matlab^® is applied to draw the fitting probability density function (PDF) curves.

Figure 9. Distribution of the TD1 outgoing velocity directions. (a) Elevation angle ϕ. (b) Azimuth angle λ of the 20 cm-sized landers. (c) Azimuth angle λ of the larger landers.

When the lander size increases, the distribution region of the elevation angle ϕ in Figure 9a moves gradually in the decreasing direction of the angle, from [59°, 72°], [37°, 42°], [30°, 37°], to [27°, 35°], respectively. The elevation angle ϕ mainly determines the trajectory’s vertical component and hopping height, simultaneously influencing the settling time and horizontal transfer distance. In Figure 9c, the lander with 30, 40, and 50 cm edge length has the angle λ distributed from [160°, 165°], [159°, 161°], to [156°, 158°] respectively. When the edge length is 20 cm in Figure 9b, the angle λ has two concentrated regions, which are [0°, 150°] and [200°, 360°], respectively. Additionally, among the four different-sized cubic landers, the angle λ of the 20 cm-sized cubic lander has the widest distribution regions, resulting in the 20 cm-sized landers’ settling positions being the most divergent. These different angles of the outgoing velocities result from the accumulation of differences in the cubic lander size during the contact.

Furthermore, the hopping times between the first landing points (TD1) and the second landing points (TD2), the settling times, and the horizontal transfer distances are depicted in Figure 10, respectively. The settling time records the hopping time from the initial release position until the lander settles, that is, the next hopping height is less than the threshold (1 cm). The horizontal transfer distance is calculated utilizing the Euclidean distances from TD1 points to the settling positions. Figure 10a shows that during the first hopping, the hopping time becomes longer from 4 to nearly 7 min when the lander edge length decreases from 50 to 20 cm. The settling time region in Figure 10b moves along the decrease direction of time with the increase of the lander size. Most 20 cm-sized landers can travel about 17 min before settling down, while the 50 cm-sized landers have the most probability of moving around 11.5 min. The elevation angle ϕ in Figure 9a leads to the changes in these settling times. When the cubic lander size increases, the elevation angle ϕ tends to be smaller. Analogous to the parabolic motion, the smaller elevation angle leads to the lower hopping height and the shorter transfer time. Similar influences on the transfer time accumulate during each collision.

Figure 10. Settling time and transfer distance distributions of different-sized cubic landers. (a) Hopping time between TD1 and TD2. (b) Settling time distribution. (c) Horizontal transfer distance distribution.

In Figure 10c, the 20 cm-sized lander with the longest transfer time has the shortest horizontal transfer distance (most distributed near 5 m). Most transfer distances of the 30 cm-sized, 40 cm-sized, and 50 cm-sized cubic landers are distributed from 10 to 15 m. When the edge length increases from 30 to 50 cm, the peak distance of PDF curve decreases from 12.08, 11.33, to 10.66 m, respectively. It indicates that the larger lander is more likely to travel shorter when the edge length increases from 30 to 50 cm. Analogous to the parabolic motion, the elevation angle ϕ is utilized to explain this relationship of the horizontal transfer distance. According to the parabolic motion, the horizontal transfer distance increases when the elevation angle of the initial velocity increases from 0° to 45° and then decreases with the angle when it is between 45° and 90°. In Figure 9a, the angles ϕ with the maximum PDF value of the 20 cm-sized, 30 cm-sized, 40 cm-sized, and 50 cm-sized cubic landers are 68.4°, 37.9°, 33.2°, and 30.8°, respectively. The differences of these angles with 45° are 23.4°, 7.1°, 11.8°, and 14.2°, respectively. The 20 cm-sized lander has the largest ϕ, which differs the most from 45°. The 30 cm-sized lander tends to have an elevation angle ϕ closest to 45°. When the lander edge increases from 30 to 50 cm, the values of angle ϕ are all less than 45°, and the difference gradually increases. Therefore, the transfer distance becomes shorter when the lander size increases. The 20 cm-sized cubic lander has the largest ϕ exceeds 45°. When ϕ is larger than 45°, the transfer distance starts to decrease theoretically. Given that the angle difference is the largest, more than 20°, the transfer distance becomes the shortest when the edge length is 20 cm.

After analyzing the trajectories when these cubic landers are released 20 m above the surfaces, the release height is increased to see its influence on the trajectory evolution. In this case, the other conditions remain the same except that the release heights increase from 20 to 30 m. The settling positions and the corresponding probability densities are illustrated in Figure 11. Similarly, 95% error ellipses of the three different-sized cubic landers are given, while only the settling positions of the 20 cm-sized landers are illustrated. When the release height increases to 30 m, the TD1 points become closer to the equator. The 20 cm-sized lander still moves along the direction of the increasing latitude, and the settling positions' latitudes become more different from the other three landers. The 95% confidence ellipses areas decrease from 0.051, 0.030, 0.026, to 0.024 when the lander edge length increases from 20 to 50 cm. The error ellipse area relationship reveals that the 20 cm lander’s trajectories are still the most divergent, while the 50 cm lander has the most gathered settling positions. Compared with the results with the 20 m release height, the areas of the error ellipses do not change significantly, indicating that the final dispersion degree is similar.

Figure 11. Final distributions of different-sized cubic landers released from 30 m.

The TD1 outgoing velocity directions, the settling times, and the transfer distances are further discussed in Figure 12. The first landing velocities’ magnitudes increase with the higher release position, while the variational trend of angle ϕ among different-sized landers in Figure 12a does not change significantly compared to Figure 9a. Larger cubic landers are still more likely to have the lower angle ϕ. Figures 12b and 12c show that the higher release position results in longer settling times and more considerable transfer distances for all-sized cubic landers, respectively. With the increase of release height, the minimum settling time of the different-sized landers increases from 10 min (Figures 10b) to 12 min, and all landers transfer larger than 5 m. However, the relationships between the settling times and the horizontal transfer distances of the different-sized landers do not change with the release position.

Figure 12. Angle ϕ, settling time, and transfer distance distributions of different cubic landers when the release height is 30 m. (a) Angle ϕ. (b) Settling time distribution. (c) Transfer distance distribution.

Comparisons about the landers with different release heights reflect that although the higher release position makes the transfer distance and time longer, it does not significantly change the settling position concentration degree. Therefore, when the pre-selected landing region is relatively larger and smoother, the spacecraft can deploy the cubic lander (such as the 40 cm-sized or 50 cm-sized) from the higher position. It may be helpful to prevent the lander from being damaged in the harsh space environment, such as the naturally ejected materials of the asteroid.

Landers With Different Shapes Cuboid Landers With Different Heights

The above cubic landers with different sizes belong to the symmetric polyhedral bodies. Their moments of inertial along the three axes of the lander frame have the same value. For analyzing the asymmetrical shape influence on the trajectory evolution, the lander shape changes from cube to cuboid in this part. Three cuboid landers are extended from the 20 cm-sized cubic lander, with one 20 cm side increased from 20 cm to 25 cm (Cuboid 1), 30 cm (Cuboid 2), and 35 cm (Cuboid 3), respectively. Masses of the landers still keep the same, and the cuboids’ masses are all uniformly distributed. The edges and the moment of inertia of Cuboid 1, Cuboid 2, and Cuboid 3 are listed in Table 2. Similar to the lander size discussions in Section 3.2, these cuboid landers are still released from the same initial position with zero velocity and no spin. Therefore, the initial release energies of these cuboid landers are kept identical. The lander shape is the only variable determining the subsequent velocity change and the trajectory evolution. Monte Carlo simulations are made for each lander shape. The 1,000 initial release states are identical to the cubic landers with the 20 m release height.

Table 2 Parameters of the Cuboid Landers

Lander	Parameters	Values
Cuboid 1	Edge length	20 × 20 × 25 cm
	Mass	10 kg
	Moment of inertia	[0.085, 0.085, 0.067]T kg·m2
Cuboid 2	Edge length	20 × 20 × 30 cm
	Mass	10 kg
	Moment of inertia	[0.108, 0.108, 0.067]T kg·m2
Cuboid 3	Edge length	20 × 20 × 35 cm
	Mass	10 kg
	Moment of inertia	[0.135, 0.135, 0.067]T kg·m2

The settling positions and the distributions of the three cuboid landers are compared with the 20 cm-sized cube in Figure 13. The dark red crosses represent the settling locations of Cuboid 1, and the settling positions of Cuboid 2, Cuboid 3, and the 20 cm-sized cube are omitted. The four 95% error ellipses are given in different colors. Like the 20 cm-sized cubic lander, the cuboid landers travel toward the higher latitude. The longitudes and latitudes of the final distributions of these cuboid landers are almost all distributed in [89.7°, 89.9°] and [4.0°, 4.3°], respectively. The corresponding 95% confidence ellipses areas are 0.028, 0.022, and 0.022 when the lander height changes from 25 to 35 cm. Comparatively, the settling points of the 20 cm-sized cubic landers are located in lower latitudes from 3.80° to 4.20° with the larger 95% confidence ellipse whose area is 0.066. It indicates that the trajectories of the three cuboid landers are more gathered than the 20 cm-sized cube. Additionally, the settling locations are more concentrated when the cuboid lander becomes more asymmetrical, that is, the lander is higher if the length and width are the same.

Figure 13. Final distributions of different-shaped cuboid landers and the 20 cm-sized cubic lander.

The comparisons of the settling times, transfer distances, and the elevation angle ϕ of these cuboid landers and the 20 cm-sized cubic lander are illustrated in Figure 14. The higher cuboid lander tends to experience the shorter time and distance during the landing motion. The settling times of the cuboid landers in Figure 14a are all <18 min. The maximum time decreases from 18, 16, to 11 min when the height increases from 25 cm (Cuboid 1) to 35 cm (Cuboid 3). Their transfer distances in Figure 14b are distributed between 5–10 m, 4–8 m, and 3–6 m, respectively, with the increase of the lander height. As for the 20 cm-sized cubic lander, it has a larger settling time between 14 and 23 min. The distribution of the transfer distance of the 20 cm cube is also broader than these cuboid landers, which can change from 1.3 to 8 m in Figure 14b. Similar to the cubic landers, the differences in the settling time and transfer distance of the three cuboid landers can be explained by the elevation angle ϕ. The elevation angle ϕ of the outgoing velocity at TD1 is illustrated in Figure 14c. The angle ϕ tends to become smaller when the cuboid height increases from 25 to 35 cm. The peak ϕ value with the largest PDF of the cuboid landers decreases from 40°, 22°, to 16°, respectively, which are all less than 45°. Therefore, the lower angle ϕ leads to a shorter settling time and smaller transfer distance in these cuboid landers. Additionally, the 20 cm cube has the largest angle ϕ, which exceeds 60°, resulting in the longest settling time.

Figure 14. Settling time, transfer distance, and angle ϕ distributions of different cuboid and cubic landers. (a) Settling time distribution. (b) Transfer distance distribution. (c) Angle ϕ.

Additionally, the PDF fitting curves of the settling time in Figure 14a and the transfer distance in Figure 14b of the asymmetrical cuboid landers (Cuboid 3, especially) have more peaks than the cubic landers. This phenomenon indicates that the asymmetrical cuboid lander tends to have more complex collision scenarios than the cubic lander. A single trajectory of Cuboid 3 is illustrated in Figure 15. This cuboid lander experiences three collisions before settling. The corresponding change of the translational energy and the rotational energy of this lander is illustrated in Figure 15b. Translational energy is much larger than rotational energy. The translational energy decreases at the collision positions, and the conversion between the two kinetic energies is observed.

Figure 15. Contact scenarios of the asymmetrical cuboid. (a) Trajectory evolution. (b) Translational energy and rotational energy of the lander.

Cylindrical Landers With Different Heights

After analyzing the cuboid lander, the corresponding cylindrical landers with identical mass and the moment of inertial are discussed in this part. The hexadecagonal-based prism (the shape of MINERVA II-1A/B of Hayabusa 2) is applied to approximate the cylindrical polyhedron model. Based on the formulas of the moment of inertial of the cylinder and the known cuboid lander inertia in Table 2, the radius and height of the cylindrical lander can be derived inversely. The radiuses of the three corresponding landers, Cylinder 1, Cylinder 2, and Cylinder 3, are solved to be 11.55 cm, with the heights are 25, 30, and 35 cm, respectively. The initial states of the Monte Carlo simulation for each cylinder remain the same. Based on the same mass and the moment of inertia, the trajectory evolution difference is exclusively influenced by the lander shape. Final distributions of these cylindrical landers are illustrated in Figure 16a. The dark red points are the settling positions of Cylinder 1. Compared with the cuboid landers, the three cylinders' settling positions overlap in an area with longitudes and latitudes are [89.7°, 90.0°] and [4.1°, 4.4°]. The 95% error ellipses areas are 0.030, 0.036, and 0.024, respectively, when the lander gradually becomes higher. It indicates that the settling positions of the three cylindrical landers are slightly divergent than the cuboid landers, and they tend to locate at higher latitudes. The elevation angles of the outgoing velocities are illustrated in Figure 16b. Compared with the cuboid lander, the angle ϕ of the cylindrical lander distributes more concentrated in a relatively smaller region from 28° to 45°. The angle ϕ with the largest PDF value of the cylindrical lander still decreases when the lander height increases. With the height increase, the most likely values of ϕ are around 42°, 35°, and 29°, respectively, larger than those of the corresponding cuboid landers (40°, 22°, and 16°).

Figure 16. Final distributions and outgoing angles of cylindrical landers with different heights. (a) Final distributions. (b) Angle ϕ.

The settling time and transfer distance of the cylindrical landers are analyzed in Figure 17. Most cylindrical landers would settle down from 10 to 16 min with the transfer distance distributed between 6 and 10 m. The time and distance distributions of the three cylindrical landers also show the characteristics of double peaks or multiple peaks because of the asymmetrical shape. Compared with the cuboid lander with the same height, the cylindrical lander tends to have a larger settling time and transfer distance due to the larger values of angle ϕ, which are closer to 45°. However, the distribution curves of these different-shaped cylindrical landers are overlapped without clear boundaries, either for the settling time or the transfer distance. Combined with the superimposed settling areas in Figure 16a, the cylindrical landers are indicated to be less sensitive to the shape change compared with the cuboid lander. This phenomenon exhibits the difference between the cylindrical surface contact and the polygon contact of the cuboid lander. It is preliminarily speculated that the sharper the angle between the side facets of the lander, the more sensitive it is to the shape change.

Figure 17. Settling time and transfer distance distributions of different cylindrical landers. (a) Settling time distribution. (b) Transfer distance distribution.

Coupled Influence of the Mass

The lander size and shape are the only variables influencing the trajectory evolution in the above parametric discussions. In this part, the masses of the 30 cm-sized, 40 cm-sized, and 50 cm-sized cubic landers are enlarged with their sizes to show the coupled influence of the lander mass. Based on the density of the 20 cm-sized homogeneous cubic lander (1,250 kg·m⁻³), the masses of the 30 cm-sized, 40 cm-sized, and 50 cm-sized cubic probes are changed from 10 to 33.75, 80, and 156.25 kg, respectively. The final distributions are illustrated and compared in Figure 18. As for one determined probe size, the settling locations of the heavier cubic probes have higher latitudes, reflecting the coupled influence of the mass.

Figure 18. Final distributions of probes with different masses. (a) A 30 cm-sized probes. (b) A 40 cm-sized probes. (c) A 50 cm-sized probes.

Figures 19 and 20 illustrate the changes of the translational energy, the rotational energy, and the total energy during the landing motion, respectively. The total energy includes the translational kinetic energy, the rotational kinetic energy, and the effective potential energy. The gravitational potential is assumed zero at the infinite distance.

Figure 19. Translational and rotational energies of the probes with the same and different masses. (a) Translational energies when the masses are identical. (b) Translational energies when the masses are different. (c) Rotational energies when the masses are identical. (d) Rotational energies when the masses are different.

Figure 20. Total energies of the probes with the same and different masses. (a) Masses are identical. (b) Masses are different.

In Figure 19, when the landers have the same masses, the change tendencies of the translational energies are identical before the first collision. The translational energy and the rotation energy interchange occurs at the first collision. Influenced by the lander size, the differences between these translational and rotational energies appear after the first collision. The two kinds of kinetic energies both decrease at the collision after the first contact. The maximum translational and rotational energies are 0.0452 and 0.0014 J, respectively. When the mass increases with the probe size, the translational energies in Figure 19b gradually become different during the descent time before the first collision. The heavier probes have the larger kinetic energy, and their landing motions last more than 15 min. By comparing Figure 19a with Figure 19c, the maximum values of the translational and rotational energy are increased to 0.75 and 0.15 J, respectively.

Adding the effective potential energy to the kinetic energy, the total energy is generated and illustrated in Figure 20, which decreases gradually with time. The total energy is negative for the lander system since the gravitational potential is assumed zero at the infinite distance. When the lander mass is identical in Figure 20a, the initial total energies are equal (about −17.42 J). The dissipation of the total energy happens at each collision. The lander position is considered unchanged during the collision, and the energy loss is mainly due to the decrease of kinetic energy. When the mass increases with size in Figure 20b, the greater the probe’s mass, the larger the absolute value of the effective potential energy at the initial time. A significant difference in the total energy exists during the entire landing motion.

In this section, the coupled influence of the lander mass on the final distribution, the kinetic energy, and the total energy is analyzed in detail. Compared with the simulations where the lander size is the only variable, the mass change significantly affects the initial energy and landing behavior. In this case, it becomes difficult to distinguish the specific influence of each factor.

Summary of Simulations

The parametric simulations are summarized in this part. The results indicate that the size and shape of a lander significantly influence its landing dynamics. Under the influence of the cubic lander size, the 40 cm-sized and 50 cm-sized cubic landers tend to have more concentrated settling locations. The transfer distance and settling time of the 50 cm-sized cubic landers are more likely to be shorter. When the release height increases from 20 to 30 m, the lander size has a similar influence on the landing trajectories despite the lander traveling longer. In terms of the different-shaped cuboid landers, the asymmetric shape makes them have more gathered settling positions than the 20 cm-sized cubic lander. Additionally, influenced by the smaller elevation angle ϕ, the more asymmetric cuboid lander is found to have more concentrated settling positions, shorter settling times, and smaller transfer distances. Furthermore, the lander shape is changed from cuboid to cylinder, with the mass and inertia identical. The cylindrical landers have the larger angle ϕ than the cuboid landers, resulting in the longer transfer distance and settling time. Additionally, the cylindrical surfaces contact makes these cylindrical landers insensitive to the shape change, compared to the polygon contact of the cuboid landers.

A lander will experience several rebounds after falling on the surfaces. During the collision period, this instantaneous process in the small body environment is very random and sensitive to the initial conditions and the lander shape/size. If the size or shape is altered, the accumulation of the collision difference over time will result in various outgoing velocity directions, establishing the dispersion degree of the subsequent trajectory evolution. Although the concrete trajectories are not constant with different simulation conditions, the outcomes in this paper can preliminarily provide some references for the lander shape and size design. In the small body missions, the lander is preferred to settle down as quickly as possible to precisely arrive at the designated area and keep communication with the mothership. Meanwhile, the lander should not be very sensitive to the initial release disturbance. Among the landers applied in this study, Cuboid 3 (the highest asymmetric cuboid lander) could settle down fastest, and the settling positions are relatively concentrated, which should be preferable. If the mission determines to utilize the symmetric cubic lander, the larger size is preferred for reducing the rebounding distance and the motion randomness when the onboard equipment mass is similar.

Conclusions

This article investigates the influence of the lander size and shape on its coupled orbit-attitude dynamical propagation. Ten different-sized landers with different shapes, classified into three main groups, are applied in numerical simulations. One thousand trajectories of the Monte Carlo simulations are implemented for the cubic landers (edge equals 20, 30, 40, and 50 cm, respectively), the cuboid landers, and the cylindrical landers (height changes from 25, 30, to 35 cm), respectively. The simulation results indicate that the size and shape of a lander significantly influence its dynamics. The larger lander size makes the cubic lander settle down faster and corresponds to more concentrated settling positions. The more asymmetric shape of the cuboid lander results in the more gathered settling positions, the shorter transfer distance, and the less settling time. When the polygon facets of the cuboid are changed into the cylindrical surfaces, the cylindrical landers settle down slower than the cuboid landers. Moreover, their settling positions become less sensitive to the lander shape change due to the cylindrical surface contact. The above results are expected to be beneficial for future lander design missions. Taking the faster settlement case as an illustration, the asymmetric cuboid lander should be first considered. When the short side is determined, the higher cuboid lander is preferable within the space scope. If the symmetric cubic lander is technically required, it is recommended to choose the larger size, leading to the shorter transfer distance.

Acknowledgments

The authors acknowledge support from the National Natural Science Foundation of China (Grant No. 11972075) and Beijing Institute of Technology (2021CX01029).

Data Availability Statement

The data sets have been uploaded to https://doi.org/10.5281/zenodo.5779077.