3D Soft-Landing Dynamic Theoretical Model of Legged Lander: Modeling and Analysis

Abstract

In this paper, a novel 3D (three-dimensional) soft-landing dynamic theoretical model of a legged lander is developed in detail as well as its numerical solution process. The six degrees of freedom motion (6-DOF) of the base model of the lander with mass center offset setting is considered in the model as well as the spatial motion (3-DOF) of each landing gear. The characteristics of the buffering force, the footpad–ground contact, and the inter-structure friction are also taken into account during the motion of each landing gear. The direct constraint violation correction is used to control the constraint stabilization of the nonlinear dynamic equation. Comparative studies between the results from the proposed model and the simulated model (built in MSC Adams) under four classical load cases show the validity of the model. Additionally, the influences of different types of contact force models, friction force models, and a friction correction model used in the soft-landing dynamic model are further investigated as a step toward understanding the soft-landing dynamic performance and the feasibility of the dynamic model method of a legged lander. The results indicate that a precise lateral force model of the footpad–ground contact is necessary to obtain the soft-landing performance of one lander during soft landing.

1. Introduction

The overall goal of China’s lunar program is to achieve China’s first manned landing on the moon by 2030 and carry out lunar scientific exploration and related technological experiments [1]. Different types of exploration equipment will be landing in the lunar polar region, which has complex lunar surface morphology and discrete lunar soil-bearing capacity characteristics [2]. Many types of modular landers with different masses, volumes, configurations, and sizes will be needed in the new-generation extraterrestrial exploration, such as landers with orientation capability [3], multifunctional landers with soft landing and locomotion [4,5,6], etc. The soft-landing dynamic theoretical model is an important method to evaluate the soft-landing performance of landers and optimize the design and arrangement of structures [7,8]. In comparison with other models (MBD (Multi-Rigid/flexible-Body Dynamic model), the FEM model, and physical test), the theoretical model has several advantages: (1) a better understanding of the load–deformation mechanism among each component, such as the footpad/ground bearing load–deformation relation; (2) easier implementation of modular designs to realize the quick design of a lander; (3) integration of multiple theoretical models to meet various design and analysis requirements; (4) reduces dependence on the commercial software such as MSC Adams or Abaqus. Therefore, it is necessary to develop a theoretical model that will satisfy the requirements of the new-generation extraterrestrial exploration and enable the easy and quick design and analysis of the next-generation lander.

In the process of soft-landing theoretical dynamic modeling, two issues require special attention: the dynamic model of the lander and the interaction force model between the footpad and the landing surface. Typically, the dynamic model of the lander comprises the main body of the lander and four landing gears (one landing gear consists of one primary and two secondary struts, an energy absorber in each strut, and one footpad) and is frequently modeled using the Newton–Euler equation method. Early works on the soft-landing dynamic theoretical model were predominantly carried out by the United States during lunar landing projects such as the Apollo mission, leading to the development of diverse theoretical models [9,10,11]. Despite the development of various theoretical models for soft-landing dynamics, many of these models were simplified, such as the simplified 3D model and two-dimensional model. As a result, the accuracy of these models may be limited in predicting the behavior of the lander during the landing process. Lavender [9] simplified the lander into a two-dimensional rigid body model with a hinged leg including elastic, damping, and crushing effects, while Alderson and Wells [10] developed a theoretical model for the Surveyor lunar lander that treated the spacecraft as a rigid body with compressible leg sets. The leg sets were treated as a plane linkage with a rigid lower link hinged to the spacecraft and a compressible, energy absorber upper link to the footpad. Zupp and Doiron [11] discussed the shortcomings of NASA’s previous soft-landing dynamics model and established a 3D dynamic simulation model. Its predictor-corrector method was the backward difference formula that maintains a specified accuracy in the integration. Takao etc. [12], based on a two-dimensional simplified landing model of a lander, proposed an improved footpad/soil-resistant theory to minimize the risk of tipping during landing. China is the third country to achieve a lunar landing, and some Chinese scholars in the field of lunar research also have conducted a large number of studies on the analytical model. Wan [13] established a 3D soft-landing model with the software MATLAB/SIMULINK and did not give the detailed calculation program flow. Yue etc. [14] built a soft-landing model of the launch vehicle using the Quasi-three-dimensional landing model and researched seven extreme landing conditions. Ke etc. [15] designed an innovative six-legged mobile lander with repetitive landing capacity and built a simplified 3D dynamic model and assessment criteria. Lin etc. [16] built a two-dimensional soft-landing theoretical model of a lunar lander, using a 7-DOF soft-landing dynamic model, and discussed the impact on soft-landing dynamic characteristics by different initial horizontal velocities, pitch angles, and inclinations of the lunar slope, etc. Yin [17] developed a planar dynamic model of a three-legged lander considering the asymmetric characteristic and the leg–leg coupling to understand the landing process of the asteroid probe with the three-legged cushioning mechanism. Yang [18] established a 3D dynamics model of China’s Mars lander considering plastic deformation parts and nonlinear contact forces. The equation of the model is too complex to be widely used in other landers. In summary, the majority of ongoing research on soft-landing dynamics were predominantly reliant on two-dimensional models, which are limited in their usage and only analyze some classical load cases. Although a few 3D models have been developed, the scalability and usability of these modeling methods are slightly insufficient, making them unable to better support future research.

The interactive model between the footpad and the lunar soil is of paramount importance, as it directly impacts the successful soft landing of the lander and its corresponding research work after soft landing. The dynamic response of the lunar soil during landings, such as the penetration depth of the footpad and the soil’s vertical and horizontal load-bearing capacity, are critical parameters in determining the stable position and angle of the lander after landing. There are three kinds of contact force models that have been derived: (1) the added mass model [19]; (2) the load model based on the Bekker theory [20], which defines contact force F as a function of the indentation depth δ: F = K × δⁿ, where K represents the stiffness parameter, and the exponent n depends on the topological properties of the contacting surfaces; (3) the simplified dynamic bearing model in the form of the function on the force versus depth and velocity [18]. Moreover, many studies of fixed-shape (non-locomoting) objects impacting and penetrating dry granular media have revealed reaction forces that can be described by

$$F=F_p\left(z\right)+\alpha v^2$$

(1)

where v and z are the velocity and depth of one object. F_p(z) is a depth-dependent force. α is the inertial drag coefficient [19]. However, during the impact, several different mechanical phenomena can occur. Tension and shear failure, localized deformation, effects of adiabatic shear, and crack propagation are only some of the important phenomena that may occur individually or simultaneously. Then, there is no uniform view on the footpad–ground bearing model when the bearing model considers other items such as penetration speed (linear and nonlinear relation), the shape of the footpad, contact type, loading weight, etc.

To meet the future design requirements and the soft-landing performance analysis of the legged lander, a 3D soft-landing dynamic theoretical model of a legged lander and its numerical solution process is developed, validated, and analyzed. In Section 2, the 3D soft-landing model of a legged lander is introduced, and the equations of kinematics and dynamics for the base model, landing gear, and footpad–ground bearing model are derived. In Section 3, the simulation program for the soft-landing model of the legged lander is developed based on the proposed method using the software MATLAB and validated by MSC Adams prototype under four classical load cases. In Section 4, different types of contact models and friction models in the footpad–ground bearing model are discussed. Furthermore, the inter-structure friction of the primary strut is also discussed.

2. Soft-Landing Model

2.1. Model Definition

A soft-landing dynamic model of the lander is shown in Figure 1. The lander consists of four landing gears, the simplified base model, and other components. Each landing gear consists of a primary strut, a footpad, and two secondary struts. Each strut is composed of outer and inner tubes connected by the sliding hinge. Moreover, the connection joints located at the points a_i, b_i, and c_i in the No. i landing gear are the universal hinge, while the joints located at the points d_i, e_i and f_i are the ball hinge.

Moreover, there are five key coordinate systems in this model: (1) the global coordinate system and local lunar coordinate system, O_g-x_gy_gz_g and O_si-x_siy_siz_si, are used to define the base reference coordinate system and the position and angle of the local lunar surface, respectively; (2) the lander coordination system, O_L-x_Ly_Lz_L, is used to define the structural distribution and motion state of the lander, where O_L is the geometry center of the lander; (3) the No. I landing gear coordinate system, O_L-x_ilgy_ilgz_ilg, is used to define the motion and load of the No. i landing gear; (4) the No. i primary strut coordinate system with origin point a_i, a_i-x_ilgy_ilgz_ilg, is used to define the position and velocity of the primary strut.

Based on the theory of spatial descriptions and transformations [21,22], the relative motion of the body or joint rotations is frequently expressed using the transition matrix in the Euler angles format. The transformation matrix ${}_i{}^{i-1}T$ consists of a translation matrix $P_i$ and a rotation matrix $R_{i-1,i}$ represented by Z-Y-X rotations, which are as follows.

$$\begin{array}{c}{}_i{}^{i-1}T={\left({}_{i-1}{}^{i}T\right)}^{-1}=\left[\begin{array}{cc}{R}_{i-1,i}& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}I& P_i\\ 0& 1\end{array}\right]=\left[\begin{array}{cc}I& P_{i-1}\\ 0& 1\end{array}\right]\left[\begin{array}{cc}{R}_{i-1,i}& 0\\ 0& 1\end{array}\right];\\ {}_i{}^{0}T={}_1{}^{0}T{}_2{}^{1}T\cdots {}_i{}^{i-1}T\end{array}$$

(2)

To express the velocity and coordinates of the landing gear conveniently, the relations of the defined coordinate systems used in the soft-landing dynamic model are defined, which is shown in Figure 2. For example, the transformation matrix ${}_{O_{LLgi}}{}^{O_L}T$ from the lander coordinate system to the No. i landing gear coordinate system can be obtained according to the rotation matrix $R_{O_L{O}_{LLgi}}$ and the translation matrix $P_{Lgi}$. The calculation equation of the rotation and translation matrix can be found in Nomenclature.

To express the soft-landing dynamic model clearly, the relation of each component in the soft-landing dynamic model is shown in Figure 3. Section 2.2 aims to calculate the position and velocity of the key point of the lander, such as the mass center of the lander M_L to use in Section 2.3, the touchdown point d_i in the No. i landing gear to use in Section 2.5, and the origin point O_Lgi of the No. i landing gear to use in Section 2.4. Moreover, Section 2.3 is used to calculate the landing gear’s velocity and length in the No. i landing gear coordinate system under the buffer and friction effect described in Section 2.6 and Section 2.7. Section 2.5 aims to calculate the binding force of the ground on the footpad, which is passed to the main body of the lander after being buffered by the landing gear.

2.2. Position and Velocity Definition of the Lander

According to the geometric relations defined in Figure 2, the coordinates of the points ${ M}_L,O_L,a_i$ under the global coordinate system, Og-x_gy_gz_g, can be denoted as follows using the calculation equations of the translation matrix in Equation (2).

$$l_M^{O_g}=\left[\begin{array}{c}{x}_m\\ y_m\\ z_m\end{array}\right];l_{O_L}^{O_g}=\left[\begin{array}{c}{x}_{O_L}\\ y_{O_L}\\ z_{O_L}\end{array}\right]; l_{a_i}^{O_g}=\left[\begin{array}{c}{x}_{ai}\\ y_{ai}\\ z_{ai}\end{array}\right];l_{d_i}^{O_g}=\left[\begin{array}{c}{x}_{di}\\ y_{di}\\ z_{di}\end{array}\right];$$

(3)

where

$$\begin{array}{c}{}_{O_L}{}^{O_g}\mathrm{T}=\left[\begin{array}{cc}{R}_{O_g{O}_L}^{}& l_{O_L}^{O_g}\\ 0& 1\end{array}\right]=\left[\begin{array}{cc}I& l_M^{O_g}\\ 0& 1\end{array}\right]\left[\begin{array}{cc}{R}_{O_g{O}_L}^{}& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}I& l_{O_L}^{M_L}\\ 0& 1\end{array}\right];\\ {}_{a_i}{}^{O_g}\mathrm{T}=\left[\begin{array}{cc}{R}_{O_g{a}_{Lgi}}^{}& l_{a_i}^{O_g}\\ 0& 1\end{array}\right]={}_{O_L}{}^{O_g}\mathrm{T}\left[\begin{array}{cc}{R}_{O_L{a}_{iLgi}}^{}& l_{a_{iLgi}}^{O_{LLgi}}\\ 0& 1\end{array}\right];{}_{d_i}{}^{O_g}\mathrm{T}={}_{a_i}{}^{O_g}\mathrm{T}\left[\begin{array}{cc}{R}_{a_{iLg}{d}_{i\mathrm{l}\mathrm{g}}}^{}& l_{d_{ilg}}^{a_{ilg}}\\ 0& 1\end{array}\right];{\mathrm{l}}_{d_{i\mathrm{l}\mathrm{g}}}^{a_{i\mathrm{l}\mathrm{g}}}=\left[d_1^i, 0, 0\right]\end{array}$$

Similarly, each point in the No. i landing gear coordinate system can be obtained. Since $l_{a_{iLgi}}^{O_{LLgi}}, l_{b_{iLgi}}^{O_{LLgi}}$ and $l_{c_{iLgi}}^{O_{LLgi}}$ are the given design parameters of the lander to define the install position of the landing gear structure, some vectors ${\mathit{l}}_{{\mathit{d}}_{\mathit{i}\mathit{L}\mathit{g}\mathit{i}}}^{{\mathit{O}}_{\mathit{L}\mathit{L}\mathit{g}\mathit{i}}},{\mathit{l}}_{{\mathit{f}}_{\mathit{i}\mathit{L}\mathit{g}\mathit{i}}}^{{\mathit{O}}_{\mathit{L}\mathit{L}\mathit{g}\mathit{i}}}$ and ${\mathit{l}}_{{\mathit{e}}_{\mathit{i}\mathit{L}\mathit{g}\mathit{i}}}^{{\mathit{O}}_{\mathit{L}\mathit{L}\mathit{g}\mathit{i}}}$ can be denoted using the generalized coordinate, $d_1^i, {\theta }_1^i, {\theta }_2^i$, according to the theory of spatial descriptions and transformations.

$$\begin{array}{l}{l}_{d_{iLgi}}^{O_{LLgi}}=f_1\left(d_1^i, {\theta }_1^i, {\theta }_2^i\right);\hspace{1em}{l}_{f_{iLgi}}^{O_{LLgi}}=f_2\left( d_{b1}^i, {\theta }_{b1}^i, {\theta }_{b2}^i\right);\\ l_{e_{iLgi}}^{O_{LLgi}}=f_3\left( d_{c1}^i, {\theta }_{c1}^i, {\theta }_{c2}^i\right);\\ l_{f_{iLgi}}^{a_{iLgi}}=f_4\left({\theta }_1^i, {\theta }_2^i\right);\hspace{1em}{l}_{e_{iLgi}}^{a_{iLgi}}=f_5\left({\theta }_1^i, {\theta }_2^i\right);\end{array}$$

(4)

According to the above equation and the principle of the virtual work, the velocity of the points $M_L,a_i,d_i,e_i,f_i$ in the global coordinate system, O_g-x_gy_gz_g, can be expressed by the Jacobian matrix and the generalized velocity, $\dot{x},\dot{y},\dot{z}, \dot{\psi }, \dot{\theta }, \dot{\varphi }, {\dot{\theta }}_1^i, {\dot{\theta }}_2^i, {\dot{d}}_1^i$, which are as follows.

$$\begin{gathered} \begin{array}{l}{v}_{M_L}^{O_g}=J_1\left(x,\mathrm{y},\mathrm{z},\mathsf{\psi },\mathsf{\theta },\mathsf{\varphi }\right)\left[\left(\dot{x},\dot{\mathrm{y}},\dot{\mathrm{z}},\dot{\mathsf{\psi }},\dot{\mathsf{\theta }},\dot{\mathsf{\varphi }}\right)\right]\\ \\ {v}_{d_i}^{O_g}=J_{di}\left(x,\mathrm{y},\mathrm{z},\mathsf{\psi },\mathsf{\theta },\mathsf{\varphi },d_1^i,{\mathsf{\theta }}_1^i,{\mathsf{\theta }}_2^i\right)\left[\left(\dot{x},\dot{\mathrm{y}},\dot{\mathrm{z}},\dot{\mathsf{\psi }},\dot{\mathsf{\theta }},\dot{\mathsf{\varphi }},{\dot{d}}_1^i,{\dot{\mathsf{\theta }}}_1^i,{\dot{\mathsf{\theta }}}_2^i\right)\right]\\ \\ {v}_{e_i}^{O_g}=J_{ei}\left(x,\mathrm{y},\mathrm{z},\mathsf{\psi },\mathsf{\theta },\mathsf{\varphi },{\mathsf{\theta }}_1^i,{\mathsf{\theta }}_2^i\right)\left[\left(\dot{x},\dot{\mathrm{y}},\dot{\mathrm{z}},\dot{\mathsf{\psi }},\dot{\mathsf{\theta }},\dot{\mathsf{\varphi }},{\dot{\mathsf{\theta }}}_1^i,{\dot{\mathsf{\theta }}}_2^i\right)\right]\\ \\ {v}_{f_i}^{O_g}=J_{ei}\left(x,\mathrm{y},\mathrm{z},\mathsf{\psi },\mathsf{\theta },\mathsf{\varphi },{\mathsf{\theta }}_1^i,{\mathsf{\theta }}_2^i\right)\left[\left(\dot{x},\dot{\mathrm{y}},\dot{\mathrm{z}},\dot{\mathsf{\psi }},\dot{\mathsf{\theta }},\dot{\mathsf{\varphi }},{\dot{\mathsf{\theta }}}_1^i,{\dot{\mathsf{\theta }}}_2^i\right)\right]\end{array} \end{gathered}$$

(5)

Moreover, the constraint equations of the geometric relationships in the landing gear can be shown in Equation (6).

$$\begin{array}{l}{l}_{f_{iLgi}}^{O_{LLgi}}=l_{a_{iLgi}}^{O_{LLgi}}+l_{f_{iLgi}}^{a_{iLgi}}=l_{b_{iLgi}}^{O_{LLgi}}+l_{f_{iLgi}}^{b_{iLgi}};\\ l_{e_{iLgi}}^{O_{LLgi}}=l_{a_{iLgi}}^{O_{LLgi}}+l_{e_{iLgi}}^{a_{iLgi}}=l_{c_{iLgi}}^{O_{LLgi}}+l_{e_{iLgi}}^{c_{iLgi}};\end{array}$$

(6)

Similarly, the velocity of point, d_iLg in the No. i landing gear coordinate system, O_lgi-x_lgiy_lgiz_lgi, can be expressed by the Jacobian matrix and the generalized coordinate velocity, ${\dot{d}}_1^i,{\dot{\theta }}_1^i, {\dot{\theta }}_2^i$, which is as follows:

$${\left[{\dot{x}}_{di}^{O_{Lgi}},{\dot{y}}_{di}^{O_{Lgi}},{\dot{z}}_{di}^{O_{Lgi}}\right]}^T=J_2{J}_3^{-1}{\left[{\dot{d}}_1^i,{\dot{d}}_{b1}^i,{\dot{d}}_{c1}^i\right]}^T=\mathit{J}{\left[{\dot{d}}_1^i,{\dot{d}}_{b1}^i,{\dot{d}}_{c1}^i\right]}^T$$

(7)

where

$$\begin{array}{c}{\left[{\dot{x}}_{di}^{O_{Lgi}},{\dot{y}}_{di}^{O_{Lgi}},{\dot{z}}_{di}^{O_{Lgi}}\right]}^T=J_2\left({\mathrm{d}}_1^i, {\mathsf{\theta }}_1^i,{\mathsf{\theta }}_2^i\right)\left[{\dot{d}}_1^i, {\dot{\theta }}_1^i,{\dot{\theta }}_2^i\right]\\ {\left[{\dot{d}}_1^i,{\dot{d}}_{b1}^i,{\dot{d}}_{c1}^i\right]}^T=J_3\left({\mathrm{d}}_1^i, {\mathsf{\theta }}_1^i,{\mathsf{\theta }}_2^i\right)\left[{\dot{d}}_1^i, {\dot{\theta }}_1^i,{\dot{\theta }}_2^i\right]\end{array}$$

Additionally, according to the principle of virtual work, the equivalent dynamic force, F_Lgi, can be obtained by the Jacobian matrix.

$$\begin{array}{c}{F}_g^T{\left[\delta x_{di}^{O_{Lgi}},\delta y_{di}^{O_{Lgi}},\delta z_{di}^{O_{Lgi}}\right]}^T=F_{Lgi}^T{\left[\delta d_1^i,\delta d_{b1}^i,\delta {\dot{d}}_{c1}^i\right]}^T;\\ F_g=J^T{Q}_{lgi}\\ F_g={\left[F_{XLgi},F_{YLgi},F_{ZLgi}\right]}^T;F_{Lgi}={\left[F_{pri}^i,F_{\mathrm{s}\mathrm{e}\mathrm{c}\_L}^i,F_{\mathrm{s}\mathrm{e}\mathrm{c}\_R}^i\right]}^T\end{array}$$

(8)

2.3. Dynamics Model of the Simplified Base Model of the Lander

Assuming that the simplified base model of the lander is a rigid body and that its total mass and the relative position of the mass center remain constant throughout the landing process, the kinematics and dynamics of the lander can be described using a 12-state variable, 6 degrees of freedom model, as outlined by the following equations [21]. Moreover, to understand the dynamic model of the whole system easily, the free-body diagram of the lander is shown in Figure 4.

$$\left\{\begin{array}{l}\dot{u}=rv-qw+F_x/m\hfill \\ \dot{v}=pw-ru+F_y/m\hfill \\ \dot{w}=qu-pv+F_z/m\hfill \\ \dot{p}=\frac{1}{I_x}\left[I_{xz}\dot{r}+I_{xy}\dot{q}+r{H}_Z-q{H}_Y+M_x\right]\hfill \\ \dot{q}=\frac{1}{I_y}\left[I_{xy}\dot{p}+I_{yz}\dot{r}+p{H}_Z-q{H}_X+M_Y\right]\hfill \\ \dot{r}=\frac{1}{I_y}\left[I_{xz}\dot{p}+I_{yz}\dot{q}+q{H}_X-p{H}_Y+M_Z\right]\hfill \end{array}\right.$$

(9)

$$\left[\begin{array}{l}\dot{x}\hfill \\ \dot{y}\hfill \\ \dot{z}\hfill \\ \dot{\psi }\hfill \\ \dot{\theta }\hfill \\ \dot{\varphi }\hfill \end{array}\right]=\left[\begin{array}{llllll}c\theta c\psi \hfill & s\theta c\psi s\varphi -c\varphi s\psi \hfill & s\theta c\psi c\varphi +s\varphi s\psi \hfill & \hfill & \hfill & \hfill \\ c\theta s\psi \hfill & s\theta s\psi s\varphi +c\varphi c\psi \hfill & s\theta s\psi s\varphi -s\varphi c\psi \hfill & \hfill & \hfill & \hfill \\ -s\theta \hfill & \mathrm{s}\varphi c\theta \hfill & c\varphi c\theta \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & 1\hfill & s\varphi \tan \theta \hfill & c\varphi \tan \theta \hfill \\ \hfill & \hfill & \hfill & 0\hfill & c\varphi \hfill & -s\varphi \hfill \\ \hfill & \hfill & \hfill & 0\hfill & s\varphi \sec \theta \hfill & c\varphi \sec \theta \hfill \end{array}\right]\left[\begin{array}{l}u\hfill \\ v\hfill \\ w\hfill \\ p\hfill \\ q\hfill \\ r\hfill \end{array}\right]$$

(10)

Typically, the forces acting on the lander include gravitational force mg, engine thrust force T, and the forces transferred from the No. i landing gear $F_{buf}^i$. These force and moment matrices, F and M, acting on the lander can be denoted as follows.

$${\mathit{F}=\mathit{R}}_{Lg}\mathit{m}\mathit{g}+\mathit{T}+{\displaystyle \sum _{i=1}^{i=4}{\mathit{R}}_{L\lg i}}{\mathit{F}}_{buffer}^i;\hspace{1em}{\mathit{M}=\mathit{R}}_{Lg}{(\mathit{l}}_{y_L}^{O_g}-{\mathit{l}}_{M_L}^{O_g})\times {\mathit{T}+\mathit{R}}_{Lg}{(\mathit{l}}_{d_i}^{O_g}-{\mathit{l}}_{M_L}^{O_g})\times {\displaystyle \sum _{i=1}^{i=4}{\mathit{R}}_{L\lg i}}{\mathit{F}}_{buf}^i;$$

(11)

where

$$\mathit{mg}=\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right];\mathit{T}=\left[\begin{array}{c}0\\ 0\\ -T\end{array}\right];{\mathit{F}}_{buf}^i={\mathit{J}}^T\left[\begin{array}{c}{N}_{pri}^i\\ N_{sce\_L}^i\\ N_{sce\_R}^i\end{array}\right]$$

$l_{y_L}^{O_g},{\mathrm{l}}_{M_L}^{O_g}$, and $l_{d_i}^{O_g}$ are the vectors expressed in the global coordinate system. The transmission force in the struts, $N_{pri}^i, N_{sec\_L}^i,N_{sec\_R}^i$, can be shown as follows:

$$N_{pri}^i=\left\{\begin{array}{cc}{F}_{pri}^i& F_{pri}^i\le F_{pri\_crush}\\ F_{\mathrm{crush}}& F_{pri}^i>F_{pri\_crush}\end{array}\right.;N_{\sec }^i=\left\{\begin{array}{cc}{F}_{\sec }^i\hfill & F_{\sec }^i\le F_{Ten\_crush\_\sec },\dot{s}\ge 0\hfill \\ F_{Ten\_crush}\hfill & F_{\sec }^i>F_{Ten\_crush\_\sec },\dot{s}\ge 0\hfill \\ F_{\sec }^i\hfill & F_{\sec }^i\le F_{\mathrm{Com}\_crush\_\sec },\dot{s}\le 0\hfill \\ F_{Com\_crush}\hfill & F_{\sec }^i>F_{Com\_crush\_\sec },\dot{s}\le 0\hfill \end{array}\right.$$

(12)

The equivalent dynamic forces $F_{pri}^i$ and $F_{sce}^i$ can be obtained by Equation (10). $F_{pri\_crush}, F_{Ten\_crush\_sec}$, and $F_{Com\_crush\_sec}$ are the crushing forces of the primary and secondary struts, which can be obtained by Equations (22)–(26). Meanwhile, the remaining driving force matrix $F_{driving}^i$ can be denoted as:

$${\mathit{F}}_{driving}^i=\left[\begin{array}{c}{F}_{driving\_pri}^{\mathrm{i}}\\ F_{driving\_\sec \_L}^{\mathrm{i}}\\ F_{driving\_\sec \_r}^{\mathrm{i}}\end{array}\right]=\left[\begin{array}{c}{F}_{pri}^i-N_{pri}^i\\ F_{\sec \_L}^i-N_{\sec \_L}^i\\ F_{\sec \_R}^i-N_{\sec \_R}^i\end{array}\right];$$

(13)

where

$$F_j^{\mathrm{i}}-N_j^{\mathrm{i}}=\left\{\begin{array}{cc}{F}_j^{\mathrm{i}}-N_j^{\mathrm{i}}& F_j^{\mathrm{i}}\ge N_j^{\mathrm{i}}\\ 0& F_j^{\mathrm{i}}<N_j^{\mathrm{i}}\end{array}\right.;j=\left[\begin{array}{ccc}pri& \sec \_L& \sec \_R\end{array}\right]$$

2.4. Dynamic Model of Landing Gear

According to the principle of the Lagrange multiplier form of dynamic equations, the dynamic equation of each landing gear can be denoted as:

$$\left[\begin{array}{cc}\mathit{M}& {\mathsf{\Phi }}_{\mathit{q}}^{\mathit{T}}\\ {\mathsf{\Phi }}_{\mathit{q}}& 0\end{array}\right]\left[\begin{array}{c} \ddot{\mathit{q}}\\ \mathsf{\lambda }\end{array}\right]=\left[\begin{array}{c}{\mathit{Q}}_{\mathit{ine}}^i+{\mathit{Q}}_{\mathit{grv}}^{\mathit{i}}\\ \mathsf{\gamma }\end{array}\right]+\left[\begin{array}{c}{\mathit{Q}}_{driving}^{\mathit{i}}\\ 0\end{array}\right]$$

(14)

${\mathit{Q}}_{driving}^{\mathit{i}}$ is the driving force in the generalized coordinate system. ${\mathit{Q}}_{\mathit{ine}}^{\mathit{i}}$ is the generalized internal force matrix, which is the velocity coupling force, and ${\mathit{Q}}_{\mathit{grv}}^{\mathit{i}}$ is the generalized gravity force matrix of the landing gear struts.

$${\mathit{Q}}_{driving}^{\mathit{i}}={\left[\begin{array}{ccccccccc}{F}_{driving\_pri}^{\mathit{i}}\hfill & 0\hfill & 0\hfill & F_{driving\_\mathbf{sec}\_L}^{\mathit{i}}\hfill & 0\hfill & 0\hfill & F_{driving\_\mathbf{sec}\_r}^{\mathit{i}}\hfill & 0\hfill & 0\hfill \end{array}\right]}^T$$

(15)

However, considering that the constraint relation Φ cannot be invalidated or deleted and is changed with the time variable during the soft landing, the Φ_qt and Φ_tt also equal zero, thus the constraint matrix γ is listed as follows.

$$\mathsf{\gamma }=-{\left({\mathsf{\Phi }}_{\mathit{q}} \dot{\mathit{q}}\right)}_{\mathit{q}}{ \dot{\mathit{q}} -2\mathsf{\Phi }}_{\mathit{qt}} \dot{\mathit{q}}-{\mathsf{\Phi }}_{\mathit{tt}}= -{\left({\mathsf{\Phi }}_{\mathit{q}} \dot{\mathit{q}}\right)}_{\mathit{q}} \dot{\mathit{q}}$$

(16)

Additionally, the numerical method for this dynamic model is direct integration with direct constraint violation correction [23,24], which can efficiently control the violations of constraint equations within any given accuracy at each time step. Compared to conventional methods such as the Newmark method or the generalized method, this method has a clear physical meaning, less calculation, an obvious correction effect, and a minor effect on the form of the dynamic equation of systems. The detailed algorithm of the method used in this paper is in Section 3.

2.5. Footpad–Ground Bearing Model

The footpad–ground bearing model comprises the vertical bearing model and the horizontal bearing model. While the vertical bearing model has been extensively studied in previous research, there is no uniform view of the calculation function, despite there being four calculation models available. To establish a widely accepted calculation model, the approximate contact calculation equation integrated into MSC Adams was utilized in this model. Meanwhile, since the surface soil exhibits low bearing capacity in the horizontal direction and is more easily movable, the horizontal bearing force is usually equal to the friction force acting on the footpad–ground interface caused by the vertical bearing force. Thus, the footpad–ground contact force vector under the local lunar coordinate system is expressed based on the geometry relation of the lander. The calculated force equation can be expressed as follows.

$$F_{lun}=\left[\begin{array}{c}{F}_X\\ F_Y\\ F_Z\end{array}\right]={\left[\begin{array}{ccc}\frac{-u{\dot{d}}_{\mathrm{locix}}}{\sqrt{{\left({\dot{d}}_{\mathrm{locix}}\right)}^2+{\left({\dot{d}}_{\mathrm{lociy}}\right)}^2}}& \frac{-u{\dot{d}}_{\mathrm{lociy}}}{\sqrt{{\left({\dot{d}}_{\mathrm{locix}}\right)}^2+{\left({\dot{d}}_{\mathrm{lociy}}\right)}^2}}& 1\end{array}\right]}^T\cdot F_Z$$

(17)

$$F_Z=\left\{\begin{array}{cc}0& d_{lociz}<0\\ K_{\mathrm{g}}{\left(d_{lociz}\right)}^n+\kappa C_g{\dot{d}}_{loci\mathrm{z}}& d_{lociz}\ge 0\end{array}\right.,\kappa =\left\{\begin{array}{cc}(d_{lociz}/d_1)& d_{lociz}<d_1\\ 1& d_{lociz}\ge d\end{array}\right.$$

(18)

2.6. Dynamic Model of the Buffers

The previous lander was equipped with multiple types of buffers in its landing gear, including the AL-honeycomb buffer, the large plastic deformation rod buffer, and the hydro-pneumatic buffer. Each buffer can be modeled as a normalized equation form, which can be expressed as follows.

$$F=a\left(\ddot{s}\right)+b\left( \dot{s}\right)+c\left(s\right)+d$$

(19)

Moreover, this equation is derived based on the principles of virtual work and the geometric relations defined by the generalized coordinates of the multi-body system. It takes the following form.

$$F_q=\frac{W\left(F,s\right)}{\partial q}=F\frac{\partial s}{\partial q}$$

(20)

Due to the slow touchdown velocity of the lander, the strain rate effect of the Al-honeycomb buffer is not considered. Then, the force of the Al-honeycomb buffer in the primary strut can be denoted as:

$$F=\left\{\begin{array}{cc}0& other\\ F_1& s_1<s\le s_0,s\le s_{his\min },\dot{s}<0\\ F_2& s_2<s\le s_1,s\le s_{his\min },\dot{s}<0\\ F_3& s\le s_2,s\le s_{his\min },\dot{s}<0\end{array}\right.$$

(21)

Similarly, the crushing force of the Al-honeycomb in the secondary strut can be denoted as:

$$F=\left\{\begin{array}{cc}0& other\\ F_{com}& s_{c1}<s,s\le s_{his\min }^s,\dot{s}<0\\ F_{Ten}& s_{t1}>s,s\ge s_{his\max }^s,\dot{s}>0\\ F_{com}^1& s_{c1}>s,s\le s_{his\min }^s,\dot{s}<0\\ F_{T\mathrm{en}}^1& s_{t1}<s,s\ge s_{his\max }^s,\dot{s}>0\end{array}\right.$$

(22)

Moreover, in the lander with the cantilever beam landing gear, the friction force between the outer and inner tubes in the primary strut must be considered. Then, the contact force can be denoted as follows:

$$\begin{array}{l}{N}_{E/F}^{\mathrm{i}}=\sqrt{{\left(F_{Ycon}^i\right)}^2+{\left(F_{Zcon}^i\right)}^2};\\ F_{con}^i=R_{{\mathit{d}}_{\mathit{i}}\mathit{iLg}}\left(\frac{{\mathit{l}}_{{\mathit{e}}_{\mathit{iLgi}}}^{{\mathit{b}}_{\mathit{iLgi}}}}{\left|{\mathit{l}}_{{\mathit{e}}_{\mathit{iLgi}}}^{{\mathit{b}}_{\mathit{iLgi}}}\right|}{\mathit{F}}_{\mathit{sec}\_\mathit{L}}^i+\frac{{\mathit{l}}_{{\mathit{f}}_{\mathit{iLgi}}}^{{\mathit{c}}_{\mathit{iLgi}}}}{\left|{\mathit{l}}_{{\mathit{f}}_{\mathit{iLgi}}}^{{\mathit{c}}_{\mathit{iLgi}}}\right|}{\mathit{F}}_{\mathit{sec}\_\mathit{R}}^i\right);\end{array}$$

(23)

The forces can be obtained by Equation (10). The generalized crush force of the primary strut in the No. i landing gear can be denoted as:

$$F_{pri}^i=\eta u_1{N}_{E/F}^{\mathrm{i}}+F_{pri}^i$$

(24)

Further, the calculation equation of η is listed in Section 2.7.

2.7. Correction Coefficient η

Due to the force-transmitting feature of the landing gear, the lateral denominational force of the primary strut was equal to the force acted by the secondary struts. Therefore, the correction coefficient η was only used in the contact force calculation between the outer and inner tubes in the primary strut. Figure 5 illustrates the force relation among the outer, inner, and secondary struts in one landing gear. Based on the structural characteristics and the contact behavior of the primary strut, two assumptions were given: (1) the deflection and angle of each cross-section at the landing gear are consistent with the deformation coordination relationship; (2) the contact pressure p upon the outer tube is distributed in the form of a cosine function by the inner tube, p = p_m × cos(θ), where p_m is the maximum pressure acting on the outer tube and with the same direction as N_E/F; (3) the contact angle θ between the outer and inner tubes is assumed to be π/2 since the inner diameter of the outer tube is approximately equal to the outer diameter of the inner tube. The force diagrams among the outer, inner, and secondary struts in one landing gear are shown in Figure 4.

Then, according to the above assumptions, the calculation equation of the correction coefficient η is denoted as follows:

$$\eta ={\eta }_1{\eta }_2$$

(25)

According to the force relations in the landing gear, the calculation equations for these correction factors are as follows. L_O, L₁, and L_L are the lengths of the different areas in one landing gear. A detailed introduction is shown in Figure 5.

$$\begin{array}{l}{\eta }_1=\frac{\left|N_A\right|+\left|N_B\right|}{N_{E/F}}=1+\frac{2L_{O1}}{\left(L_I+L_O-L_L\right)};\\ {\eta }_2=\frac{{\displaystyle \int \mathit{p}dl}}{N_{E/F}}=\frac{{\displaystyle {\int }_0^{\pi /2}{p}_m{l}_R\cos {\theta }_c}d{\theta }_c}{{\displaystyle {\int }_0^{\pi /2}{p}_m{l}_R\cos {\theta }_c}\cos {\theta }_{c}d{\theta }_c}=\frac{4}{\pi }\end{array}$$

(26)

3. Simulation and Model Validation

3.1. Program Flow

A four-legged lunar lander with a total mass of 16 tons was established as a demonstrative application to validate the dynamic model used in this paper, based on our previous research [25]. The numerical simulation was performed using MATLAB R2021a on a desktop computer with a 5.3 GHz CPU and 64 GB RAM. With a time step of 0.000005 s for a total simulation setting of 1 s, so the model only required 290 s to complete. Moreover, the numerical method in this dynamic model was the direct integration with the direct constraint violation correction to control the constraint stabilization. The calculation steps in this model code can be shown as follows.

Step 1: at time t₀ = 0, define the initial values of the designed variable:

The design variables of the lunar lander include the initial position and velocity matrix, mass property values of the lander, design parameters of the landing gear such as the crushing force of each strut, and the install position and length parameters of one landing gear. Additionally, the loading coefficient of the soil, the friction coefficient between the structure and ground or structure, as well as the value of the time increment and total simulation time, are also considered design variables.

Step 2: at time t_n, identify whether the footpad of the No. i landing gear touched down on the relative local surface or not:

Based on the generalized position and velocity of the No. i landing gear, and the position and velocity vector (at time t_n) of the lander, calculate the position and velocity matrix of point d_i, ${\dot{d}}_i$ in the No. i local lunar surface coordinate system using the Jacobian matrix in Equation (5) and the translation matrix in Equation (2). If d_lociz < 0, the transmission force of the No. i landing gear is equal to zero, then perform step 5; if touch down happens, obtain the footpad–ground bearing force vector ${\mathit{F}}_{\mathit{l}\mathit{u}\mathit{n}}^{\mathit{i}}$ using Equations (17) and (18). Based on the Jacobian matrix in Equations (7) and (8), calculate the equivalent dynamic force matrix ${ F}_{\mathrm{lgi}}$ in the No. i landing gear and the component force in each strut.

Step 3: at time t_n, identify whether the buffer of the No. i landing gear is crushed or not:

According to the current value of the length and velocity in each strut at time t_n, obtain the crushing force of the Al-honeycomb buffer of each strut using Equations (21)–(26). According to Equations (12)–(13), the transmission force in the struts and the driving force of the relative struts can be obtained. Using the driving force of the relative strut, the dynamic model of one landing gear is calculated using Equations (14)–(16). In addition, using the transmission force of the relative strut, the forces and moments acting on the lander can be obtained from the No. i landing gear using Equation (11). When the equivalent dynamic force is no less than the defined crush force of the Al-honeycomb buffer of the strut, the buffer begins to crush. The transmission force in the strut is equal to the crushing force of the relative strut. Moreover, the driving force in the relative strut is equal to the remaining force. However, if the equivalent force is less than the crushing force, the buffer does not crush. The driving force in the relative strut is equal to zero, and the transmission force in the strut in the No. i landing gear is equal to the equivalent dynamic force.

Step4: at time t_n, calculate the dynamic model of the landing gear i and check the constraint stabilization:

Solve the variables in Equation (11), then the position and velocity ${\widehat{q}}^{\mathrm{n}+1}$, ${\widehat{\dot{q}}}^{\mathrm{n}+1}$ are obtained in turn using the direct integration method. Check the constraint stabilization of the nonlinear equation, if the result satisfies the constraint stabilization equation: $0<‖{\Phi }_i‖<10^{-10}$, $0<‖{\Phi }_{q}q‖<10^{-10}$, then $q^{\mathrm{n}+1}={\widehat{q}}^{\mathrm{n}+1},{\dot{\mathrm{q}}}^{\mathrm{n}+1}={\widehat{\dot{q}}}^{\mathrm{n}+1}$, go to step 5.

Otherwise, firstly, calculate $\mathsf{\Delta }\mathit{q}=-{\Phi }_{\mathit{q}}^{+}\Phi ;$ ${\Phi }_{\mathit{q}}^{+}={{\Phi }_{\mathit{q}}}^T{\left({\Phi }_{\mathit{q}}{{\Phi }_{\mathit{q}}}^T\right)}^{-1}$ to carry out the displacement constraint violation correction, then obtain $q^{\mathrm{n}+1}={\widehat{q}}^{\mathrm{n}+1}+\mathsf{\Delta }q$, check $0<‖{\Phi }_i‖<10^{-10}$ again, and, if not satisfied, perform the above calculation process again until the equation satisfies. Secondly, using the $q^{\mathrm{n}+1}$ obtained from the above displacement correction, perform the velocity constraint violation correction $\mathsf{\Delta }\dot{q}=-{\Phi }_{\mathit{q}}^{+}{\Phi }_q\widehat{\dot{q}}$ then obtain ${\dot{q}}^{\mathrm{n}+1}={\widehat{\dot{q}}}^{\mathrm{n}+1}+\mathsf{\Delta }\dot{q}$, check $0<‖{\Phi }_{q}q‖<10^{-10}$ again, and, if not satisfied, perform the above calculation process again until the equation is satisfied. At last, the position $q^{\mathrm{n}+1}$ and ${\dot{q}}^{\mathrm{n}+1}$ after constraint violation correction will be taken in the next dynamic model.

Step 5: identify whether all landing gears have been calculated or not:

If the transmission force of all the landing gears has been calculated, then go to step 6. Otherwise, return to step 2.

Step 6: at time t_n, calculate the dynamic model of the lander:

Based on the forces and moments from all the landing gear systems, solve the acceleration variable of the lander using Equations (9) and (10) and obtain the position and velocity vector (at time t_n+1) of the lander in the lander coordinate system in turn using the direct integration. Then, based on the transient matrix theory, calculate the position and velocity value of the lander in the global coordinate system.

Step 7: repeat steps 2 to 7 until the total time is over:

To clearly introduce the detailed calculation steps for the soft-landing model, the flow chart of the touch down dynamic model of the lander during a soft landing is presented in Figure 6.

3.2. Program Flow

Figure 7 and

Table 1. The parameter information of a lunar lander.

Mass Property		C.L. of Points		Crush Force of Buffers				Contact Property
				Pri. Strut		Sec. Strut
Ix	8.18 × 104	$a_1$	(3.89, 0, −2.18)	F1	78,000	Fcom	43,000	Kg	(1000)1.5 ×105
Iy	8.18 × 104	$b_1$	(3.37, 1.23, −1.01)	F2	156,000	Ften	30,000	Cg	10,000
Iz	8.17 × 104	$c_1$	(3.37, −1.23, −1.01)	$s_1$	0.2	$s_{t_1}$	0.27	$D_1$	0.0001
Mass	1.6 × 104	$d_1$	(5.48, 0, 0.94)	$s_2$	0.475	$s_{c_1}$	0.27	us	0.4/0.4
		$e_1$	(4.65, 1.09, −0.505)					$u_{\mathit{d}}$	0.1
		$f_1$	(4.65, −1.09, 0.505)					$u_1$	0.1
		MC	(0, 0, −3.818)					η	0
								vs	0.1
								vd	1

Note: (1) unit: mass: kg, length: m, velocity: m/s, angle: deg, force: N, moments of inertia: kg×m²; (2) Kg, Cg, $D_1$, and u_s, $u_{\mathit{d}}$, v_s, $v_{\mathit{d}}$ are the default values in MSC Adams.

present the configuration and parameters of a legged lander. Further details on the parameters and modeling of the lander can be found in [25]. To analyze the soft-landing dynamics of a legged lander, four severe load cases [25] are chosen in

Table 2. The load case information for soft-landing analysis.

Load Cases	Velocity m/s		Attitude Angles /Deg	Lunar Surface Slope Angle /Deg
	Vertical	Horizontal
LC-1	4	−0	0/0/0	0
LC-2	4	1	45/0/0	0
LC-3	4	1	45/−4/0	8
LC-4	4	1	0/4/0	8

Note: (1) unit: mass: kg, length: m, velocity: m/s, angle: deg. (2) Attitude angles are the rotation angle relation to the global coordinate system in terms of Z-Y-X.

, including vertical loading load case (LC-1), the high-overloading load case (LC-2), the easily turnover load case (LC-3), and the maximum compression of primary strut load case (LC-4). To conduct a comparative analysis, two analytical models for the lander were generated using the methodology proposed in this paper and MSC Adams. Notably, there are three main differences between the two models, consisting of (1) the contact and friction force calculation equation differs: the STEP functions used in Adams are non-linear and non-clearly expression equations, while the functions used in the proposed method are linear functions; (2) the defined method of the crushing force: in MSC Adams, the widely used method of the SFROCE function was used, only considering the velocity’s direction and the crushing displacement of the buffer, without the maximum/minimum history displacement; (3) the numerical calculation method: the Runge–Kutta numerical method was used in MSC Adams, while the direct integration with the direct constraint violation correction was used in the proposed model.

Since the crushing length of each buffer after landing is the most direct and effective way to validate the simulated results [18], the crushing length of each strut after touch down is thereby analyzed emphatically in this work as well as the deviation data of the crushing length.

Table 3. The crush length of each Al-honeycomb buffer during soft landing.

No.	LC-1			LC-2			LC-3			LC-4
	Theory	Adams	Dev	Theory	Adams	Dev	Theory	Adams	Dev	Theory	Adams	Dev
1 1	287.0	288.8	−1.8	314.3	319	−4.7	360.3	354.1	9.8	390.0	385.0	5.0
2 2	49.8	50.3	−0.5	3.9	5.1	−1.2	−112.1	−132.2	19.1	4.4	5.0	−0.6
3 3	49.8	50.7	−0.5	72.6	74.5	−1.9	77.0	74.2	3.8	4.4	5.1	−0.7
4 4	287.0	288.7	−0.7	244.0	240.0	4.0	255.3	252.7	2.6	280.6	261.0	19.6
5 5	49.8	50.0	−0.2	41.2	43.5	−1.8	46.1	49.4	−3.3	7.3	7.3	0
6 6	49.8	50.3	−0.5	113.7	119.7	6.0	132.3	135.6	3.3	151.9	164.4	−13.5
7 7	287	289.2	2.0	244.0	242.9	1.1	255.4	251.2	4.2	277.9	261.0	16.9
8 8	49.8	50.6	−0.8	113.7	121.4	−7.7	132.6	135.6	3.0	98.3	121.9	−20.9
9 9	49.8	50.7	−0.9	41.2	45.5	−4.3	46.0	49.4	3.0	98.3	121.5	−20.9
10 10	287.0	287.5	2.0	314.2	319.5	−5.3	360.3	354.0	6.3	280.6	261.6	19.0
11 11	49.8	51.7	−1.9	72.6	75.1	−2.5	77.0	74.5	2.5	151.9	165.0	−14.1
12 12	49.8	51.7	−1.9	3.9	5.2	−1.3	−114.7	−132.2	18.5	7.3	7.3	0

¹ No. 1 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 1 landing gear. ² No. 2 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 1 landing gear. ³ No. 3 measure option: the crush length of each Al-honeycomb in the secondary strut (right) in the No. 1 landing gear. ⁴ No. 4 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 2 landing gear. ⁵ No. 5 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 2 landing gear. ⁶ No. 6 measure option: the crush length of each Al-honeycomb in the secondary strut (right) in the No. 2 landing gear. ⁷ No. 7 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 3 landing gear. ⁸ No. 8 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 3 landing gear. ⁹ No. 9 measure option: the crush length of each Al-honeycomb in the secondary strut (right) in the No. 3 landing gear. ¹⁰ No. 10 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 4 landing gear. ¹¹ No. 11 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 4 landing gear. ¹² No. 12 measure option: the crush length of each Al-honeycomb in the secondary strut (right) in the No. 4 landing gear.

and Figure 8 show the crushing length of each strut in the four classical load cases as well as their error data under the theoretical dynamic model and multi-body model in MSC Adams. The comparative results show that the crushing length of each strut from the theoretical dynamic model was close to that from MSC Adams in general. Most of the error ratios of the results ranged between 2.0% and 7.0%. The max deviation crushing length was nearly 20.9 mm at the primary strut of LG-3 (No. 3 landing gear) in LC-4, and the relative error ratio was 17.0% based on the crushing length from MSC Adams that was 261.6 mm. The main reason for this calculation deviation is the differences between the calculation equations for contact and friction force, especially the calculation of friction force. Moreover, the other reason may be the SFROCE function used in the definition of the buffer force between the software and the theoretical model, which causes a small error during landing. The SFROCE function defined in MSC Adams does not consider the maximum history displacement of each strut. Once the reciprocating motion of a strut occurred during landing, some errors were taken in.

Furthermore, it also can be found that the max error ratios were 23% at LC-2 based on the fact that the crushing length from MSC Adams was 5.2 mm. Its relative deviation was 1.7 mm at the secondary struts in LG-1 and LG-4. The reason may be the different numerical solution process between the software and the theoretical model, which causes a small error during landing. In the numerical method, as we all know, there is, more or less, some deviation between the results from each other method. Moreover, according to Figure 9, it is easily found that the stroke of each strut in the lander from two models have good agreement under four load cases. In conclusion, the agreement of the results proves the availability of soft-landing prediction of the theoretical dynamic model proposed in this paper.

4. Discussions

4.1. Different Kinds of the Footpad–Ground Bearing Models

Although there are different types of footpad–ground bearing models, there is no uniform view of the functional form of the footpad–ground bearing model. The essence of the footpad–ground bearing model is a dissipative contact force model, which is a pivotal tool to predict the contact force and energy dissipation characteristics of the soil during soft landing. The soft-landing performance of a lander is therefore highly dependent on the accuracy and precision of the contact force model. To better understand the differences and similarities between contact force models, some frequently used impact contact force models, such as the Hertz contact model [26], Hertz contact + linear damping factor model [27], Hertz contact + step damping factor model [28,29] Hertz contact + bilinear damping factor model [19], and Hertz contact + hysteresis damping factor model [30], were discussed in this section, which are listed in

Table 4. Some commonly used contact force models defined in this work.

Contact Model Types	Equation	Contact Model	Equation
Hertz contact	$F=K{\delta }^n$	Hertz contact + linear damping	$F=K{\delta }^n+D\dot{\delta }$
Hertz contact + step damping	$F=\left\{\begin{array}{cc}K{\delta }^n+D\dot{\delta }& \delta \le d\\ K{\delta }^n& \delta \ge d\end{array}\right.$	Hertz contact + bilinear damping	$F=K{\delta }^n+D_1\dot{\delta }+D_2{\dot{\delta }}^2$
Hertz contact + hysteresis damping 1	$F=K{\delta }^n+\chi {\delta }^{{}^n}\dot{\delta }$

¹ $\chi =\frac{3\left(1-c_r\right)}{2{\dot{\delta }}^{\left(-\right)}}K,c_r=1-\alpha {\dot{\delta }}^{\left(-\right)}.$ $\alpha$ is a constant value ranging between 0.008 and 0.32.

. According to Table 4, it was easy to see that the difference between the five contact models was the damping option, aimed to accurately describe energy dissipation in the collision process. Detailed advantages and disadvantages for most contact models are discussed in reference [31].

Table 5. The parameter information of the five different contact force calculation models.

Model Name 1	Model Type	Parameter Defined in the Model
Hertz contact	Hertz contact	K = (1000) 1.5 × 1.0 × 105, n = 1.5
Kelvin–Voigt model	Hertz contact + linear damping factor	K = (1000)1.5 × 1.0 × 105, D = 10 × 1000, n = 1.5
Kelvin–Voigt 1	Hertz contact + step damping factor	K = (1000)1.5 × 1.0 × 105, D = 10 × 1000, n = 1.5
Kelvin–Voigt 2	Hertz contact + bilinear damping factor	K = (1000)1.5 × 1.0 × 105, D1 = 10 × 1000, D2 = 10 × 1000, n = 1.5
Hunt–Crossley	Hertz contact + hysteresis damping factor	K = (1000)1.5 × 1.0 × 105, n = 1.5

¹ The hertz contact + step damping factor model and the Hertz contact + bilinear damping factor model are the contact models with relation to velocity options, which are similar to the Kelvin–Voigt model. To clearly express the difference between the contact models, the names Kelvin–Voigt 1 and Kelvin–Voigt 2 are defined as the two contact models, respectively.

shows the parameter information of the five different contact force calculation models. To discuss the effect of the five different contact force models on the soft-landing performance, the parameters of each contact force model were chosen based on the rigid surface contact theory of MSC Adams. Moreover, the friction calculation method was the penalty function method, and the friction coefficient u_s = u_d = 0.4.

Table 6. The crush length of each Al-honeycomb buffer under five contact models.

Info.	No.	Base Result	Hertz	Kelvin–Voigt	Kelvin–Voigt 1	Kelvin–Voigt 2	Hunt–Crossley
LC-1	1 1	287.0	288.8	287.3	288.8	287.0	289.8
	2 2	49.8	50.5	49.8	50.5	49.8	50.7
	3 3	49.8	50.5	49.8	50.5	49.8	50.7
LC-2	4 4	314.0	312.9	314.1	311.4	314.0	317.7
	5 5	3.9	6.3	4.0	6.4	3.9	12.1
	6 6	41.1	41.6	41.2	41.8	41.1	41.6
	7 7	244.0	244.5	244.1	244.9	244.0	244.1
	8 8	113.6	112.7	113.9	114.0	113.6	109.1
	9 9	73.4	74.2	72.7	75.1	73.4	74.0

¹ No. 1 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 1 landing gear. ² No. 2 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 1 landing gear. ³ No. 3 measure option: the crush length of each Al-honeycomb in the secondary strut (right) in the No. 1 landing gear. ⁴ No. 4 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 1/4 landing gear. ⁵ No. 5 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 1 landing gear. ⁶ No. 6 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 2 landing gear. ⁷ No. 7 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 2/3 landing gear. ⁸ No. 8 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 3 landing gear. ⁹ No. 9 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 4 landing gear.

and Figure 10 present the error data and crushing length of each strut under five different contact models for two classical load cases (LC-1 and LC-2). The simulation results, based on a comparative analysis of each buffer in the lunar lander after landing, indicated that the crushing length of each strut in the theoretical dynamic model under different contact force models was generally similar. The deviation ratios of most results under the five contact force models were less than 5%. The maximum deviation ratio of crushing length was nearly 200% at the secondary struts of L.G.-1 in LC-2 under the H–C contact force model, and the crushing value was 12.1 mm. The reason is due to the differences in the contact force calculation equations. To better understand these differences, the results of LC-1 were analyzed, primarily focusing on the deviation of the contact force time history curves of the footpad–ground bearing model, as shown in Figure 11. The results in Figure 11 indicated that the Hunt–Crossley contact model, the Hertz contact model, and the Kelvin–Voigt 1 model were the most influential models on the footpad–ground bearing force, with many oscillations in the load time history curves of these three models. This is the main reason for the deviation in the crushing length results of each strut. In contrast, the load time history curves of the Kelvin–Voigt and Kelvin–Voigt 2 models have fewer oscillations. Compared to the results of the Hertz contact model, the reason for this difference may be the effect of the velocity option in the calculation equation of the contact force model, which is usually used as a damping option to describe the energy dissipation and suppress data fluctuations during the collision process.

4.2. Friction Analysis

Friction is important in various fields such as engineering, physics, and materials science. Usually, the magnitude of the friction force depends on the normal force pressing the two surfaces together and the coefficient of friction between the two surfaces. The friction coefficient is a dimensionless constant that represents the friction characteristics of the two surfaces. It depends on various factors such as the nature of the two surfaces in contact, the roughness of the surfaces, the temperature, the relative speed of the surfaces, etc. To understand more clearly the lateral interaction of the footpad–ground bearing model, the default contact force model shown in Equations (17) and (18) were chosen first. Meanwhile, six commonly used friction models [32,33] were discussed in this section. The force calculation equations for the six classical friction models used in the engineering field are shown in

Table 7. The force calculation equations of the six different friction calculation models.

Friction Model Types 1	Calculation Equations 2	The Change Rule Figures
The static Coulomb friction model	$F_C=\left\{\begin{array}{cc}{u}_{s}N\mathrm{sgn}(v)& v\ne 0\\ \min (F_e,u_{s}N)\mathrm{sgn}(v)& v=0\end{array}\right.$
The static Coulomb friction model with viscous effect	$F_C=\left\{\begin{array}{cc}{u}_{s}N\mathrm{sgn}(v)+C_{v}v& v\ne 0\\ \min (F_e,u_s^{}N)\mathrm{sgn}(v)& v=0\end{array}\right.$
The regularized Coulomb friction model	$F_C=\left\{\begin{array}{cc}{u}_{s}N\left(\frac{1-e^{-\frac{3\left|v\right|}{v_0}}}{1-e^{-3}}\right)\mathrm{sgn}(v)& \left|v\right|\le v_0\\ u_{s}N\mathrm{sgn}(v)& \left|v\right|>v_0\end{array}\right.$
The stiction + dynamic Coulomb model	$F_C=\left\{\begin{array}{cc}{u}_{s}N\mathrm{sgn}(v)& v\ne 0\\ \min (F_e,u_{d}N)\mathrm{sgn}(v)& v=0\end{array}\right.$
The stiction + Stribeck + Coulomb + viscous friction model	$F_f=\left(F_D^{}+\left(F_S^{}-F_D^{}\right)e^{-{\left(\frac{\left|v\right|}{v_s}\right)}^{\delta }}\right)\mathrm{sgn}(v)$
The stiction + modified Stribeck + Coulomb + viscous friction model	$F_f=\left\{\begin{array}{cc}\left(-\frac{F_S}{v_0^2}{\left(\left|v\right|-v_0\right)}^2+F_S\right)\mathrm{sgn}(v)& \left|v\right|<v_0\\ \left(F_C^{}+\left(F_S^{}-F_C^{}\right)e^{-\epsilon {\left(\frac{\left|v\right|}{v_s}\right)}^{\delta }}\right)\mathrm{sgn}(v)& \left|v\right|\ge v_0\end{array}\right.$

¹ Stribeck effect is the effect that describes the decrease in the friction force as the relative tangential velocity increases; the viscous effect can be established as a proportion of the relative tangential velocity of the sliding surfaces, F = u × v; the stiction effect is the manifestation that static friction is greater than dynamic friction mechanics; ² v_o is the stiction velocity, and when v is no greater than v_o, the static force is the main force in friction force.

.

Table 8. The parameter information of the six different friction calculation models.

No.	Friction Model	Parameter
TY-1	Static Coulomb friction	$u_s$ = 0.4
TY-2	Static Coulomb friction model with viscous effect	$u_s$ = 0.4; $C_v$ = 100
TY-3	Regularized Coulomb friction	$v_0$ = 0.1 m/s; $u_s$ = 0.4
TY-4	Stiction + dynamic Coulomb	$u_d$ = 0.4
TY-5	Stiction + Stribeck + Coulomb + viscous friction	$v_s$ = 1 m/s; $\delta =0.85$
TY-6	Stiction + modified Stribeck + Coulomb + viscous friction	$v_0=v_s=\frac{0.1m}{s}; \epsilon =0.85$ $\delta =0.85$

provides the parameter information for the six different friction calculation models.

Table 9. The crush length of each Al-honeycomb buffer after soft landing.

	Measure Option	Base Result	TY-1	TY-2	TY-3	TY-4	TY-5	TY-6
LC-1	1 1	287.0	295.1	295.0	290.6	295.1	295.1	290.2
	2 2	49.8	48.4	48.4	49.3	48.4	48.4	49.5
	3 3	49.8	48.4	48.4	49.3	48.4	48.4	49.5
LC-2	4 4	314.0	368.3	324.4	316.0	368.3	323.7	313.4
	5 5	4.0	10.2	7.4	7.2	10.2	7.4	6.8
	6 6	41.1	39.6	39.8	40.9	39.6	39.8	41.2
	7 7	244.0	252.7	249.0	245.3	252.7	249.1	245.5
	8 8	113.6	92.4	106.0	109.9	92.4	107.0	111.8
	9 9	73.4	62.0	71.8	75.7	62.0	73.2	76.4

¹ No. 1 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 1 landing gear. ² No. 2 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 1 landing gear. ³ No. 3 measure option: the crush length of each Al-honeycomb in the secondary strut (right) in the No. 1 landing gear. ⁴ No. 4 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 1/4 landing gear. ⁵ No. 5 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 1 landing gear. ⁶ No. 6 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 2 landing gear. ⁷ No. 7 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 2/3 landing gear. ⁸ No. 8 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 3 landing gear. ⁹ No. 9 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 4 landing gear.

and Figure 12 show the deviation ratio of the crushing length and the crushing length of each strut under six different friction models in the two classical load cases. Based on the comparative results of each strut in the lander, it is shown that the crushing length of each strut from the theoretical dynamic model under different friction force models was close to each other. The max deviation for crushing length was 54.3 mm at LC-2 under TY-1 and TY-4 models, and the relative error was 17.3% at the primary struts, while the relative value of the crushing length from the theory model was 312.4 mm. The max deviation ratio for crushing length was 18.7% at LC-2 under TY-1 and TY-4 models, and the relative deviation was 21.2 mm at the secondary struts, while the relative value of the crushing length from the theory model was 113.6 mm. Causing these calculation deviations was the velocity option defined in the friction calculation model and its coupling cumulative deviation mechanism. Since the TY-1 and TY-4 friction models do not consider the constraint of the speed option, its component force in the direction of the main strut is less than that of the other friction force models. Meanwhile, due to less influence on the friction model by the lateral velocity of the lander in LC-1, the results among each friction model were relatively closer than those in LC-2. In addition, compared with the results of the contact force model, it was found that the influence of the friction model on the results was greater than that of the contact force model. Therefore, more attention should be paid to the friction force model in the footpad–ground force model during the modeling and analysis process of the soft-landing dynamics of the lander.

4.3. Correction Coefficient η

According to the force-transmitting feature of the landing gear, the lateral denominational force of the primary strut equals the force acted by the secondary strut. The correction coefficient η is thereby only used in the contact force calculation between the outer and inner tube in the primary strut. To analyze the effect of the correction factor η, the comparison analysis of the crushing length and its deviations were discussed in this section.

Table 10. The crush length of each Al-honeycomb buffer during soft landing for correction factor.

Info.	Measure Options	$\mathbf{Without} \mathit{\eta }$/mm	$\mathbf{Within} \mathit{\eta }$/mm	Error/%
LC-1	1 1	287.0	276.9	3.6
	2 2	49.8	48.2	3.3
	3 3	49.8	48.2	3.3
LC-2	4 4	314.0	300.9	4.4
	5 5	4.0	6.6	−39.0
	6 6	41.1	40.0	2.7
	7 7	243.7	235.2	3.6
	8 8	113.6	109.9	3.4
	9 9	73.4	73.3	0.1

¹ No. 1 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 1 landing gear. ² No. 2 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 1 landing gear. ³ No. 3 measure option: the crush length of each Al-honeycomb in the secondary strut (right) in the No. 1 landing gear. ⁴ No. 4 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 1/4 landing gear. ⁵ No. 5 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 1 landing gear. ⁶ No. 6 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 2 landing gear. ⁷ No. 7 measure option: the crush length of each Al-honeycomb in the primary strut in the No. 2/3 landing gear. ⁸ No. 8 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 3 landing gear. ⁹ No. 9 measure option: the crush length of each Al-honeycomb in the secondary strut (left) in the No. 4 landing gear.

and Figure 13 show the deviation data of the crushing length and the crushing length of each strut with and without correction factor η in the two classical load cases. According to the comparative results of each buffer, the simulation results showed that the crushing length of each strut in the theory dynamic model with and without friction correction was close to each other in general. The max error in crushing length was 13.1 mm at the primary strut 1/4 in LC-2, and the relative error ratio was 4.4%. Causing this calculation deviation was mainly the added friction in the crushing force of the primary strut, which is mainly led by the lateral force. In fact, the friction force upon the primary strut is decided by two factors, namely the friction coefficient and the lateral force in the landing gear. Since the friction coefficient is only dependent on the material properties, the added friction is thereby only decided by the lateral force in the landing gear. According to Equations (23) and (24), the lateral force upon the primary strut was mainly decided by the defined crushing force and the relative velocity of the secondary strut.

The deviation ratios of the results in the secondary strut ranged between 0.1% and 4.4%. However, there is a large deviation ratio among the results, 39% at the primary strut in LC-2. Considering that the base value is 4.0 and some deviation may exist, the result can be accepted. In conclusion, the agreement of the results proved the friction force upon the primary strut has few effects on the soft-landing performance of the legged lander.

5. Conclusions

A novel 3D soft-landing dynamic theoretical model of a legged lander is developed in detail as well as its numerical solution process. The six degrees of freedom motion (6-DOF) of the base model of the lander with mass center offset setting is considered in the model as well as the spatial motion (3-DOF) of each landing gear. The characteristics of the buffering force, the footpad–ground contact, and the inter-structure friction are also taken into account during the motion of each landing gear. Some salient conclusions are summarized as follows:

Comparative analysis between the theory dynamic model and multi-body model in MSC Adams under four classical load cases were carried out. The results show that the crushing length of each strut from the theory dynamic model is close to that from the MSC Adams Software in general. Despite there being some max deviation error in some struts, the theoretical dynamic model established in this paper remains feasible due to the differences in the contact and friction force calculation models as well as the solving method in each numerical solution process.

Comparative analysis among the five classical contact models and the six commonly used friction models shows that the crushing lengths of most struts from the theory dynamic model are relatively close to each other. However, it is in the friction model result that the max deviation of crushing length is 53.8 mm at the primary strut in LG-3 at LC-2, and the average relative error is 17.3%. The relative value from the theoretical model is 312.4 mm. Causing the deviation error is the difference between the velocity option in the friction models. Therefore, building a precise friction model in the footpad–ground bearing model during the soft-landing process is necessary to obtain the soft-landing performance of one lander.

The friction correction between the outer and inner tubes in the primary strut is also discussed. The results show that the deviation of the model with and without the friction correction coefficient is not significantly obvious.