Numerical solution, simulation and testing of the thermal dynamic characteristics of ball-screws

Abstract

This research focuses on the characteristics of how a ball-screw changes its temperature corresponding to the periodic change of the end-most heat source. Based on the theory of heat transfer, this paper supplies a numerical solution to the non-homogeneous equation of heat transfer through the group explicit (GE) finite difference approach. The temperature distribution for the ball-screw at different times and directions is shown. By simulating and modeling the temperature field and thermal deformation of the ball-screw under periodically varying heat sources, this paper describes the thermal dynamic characteristics of the ball-screw under such conditions. By testing the thermal dynamic characteristics of the ball-screw, the numerical solution is validated.

1 1 Introduction

The probers, including lander and rover, play a pivotal role in scientific exploration on extraterrestrial bodies, such as Moon and Mars [1]. Luna-2, the first hard lander to reach lunar surface, can only perform a transient probe before colliding with the moon. Probers with landing capability are produced to realize surface exploration on the Moon; they usually adopt three or four foldable legs with invert-tripod mechanism, such as Surveyor-1 [2] and Luna-16 [3], or with cantilever mechanism, such as Apollo project [4] and Chang’e 3, 4, and 5 [5–7]. For probers with invert-tripod or cantilever mechanism, their leg often consists of a main branch chain and two multifunctional auxiliary branch chains. All branch chains use internal aluminum honeycomb material to absorb landing impact energy. However, the exploration capacity of such probers is still limited by their immovable nature. The drill or dig devices [8] mounted on landers can only sample the soil surrounding the fixed landing site, which may possibly be polluted by engine propellant. Therefore, researchers have proposed two main strategies to improve the exploration capacity of probers.

The first strategy is to use the cooperative work mode of lander and rover, such as the cooperation of Apollo lander and lunar roving vehicle (LRV) [9], and Chang’e lander and Yutu-1 and Yutu-2 [5,6]. However, all the existing rovers, such as LRV and Yutu-1 and Yutu-2, are wheeled mobile robots. Although they possess excellent mobility on relatively even terrain, their traversing performance on irregular terrains, such as gully, slope, and rock, remains poor even under active suspension system. At present, the legged robots, which are driven by hydraulic actuator or electric motor, attract increasing attentions owing to the excellent mobility performance on many challenging terrains [10]. Compared with the hydraulic actuator that is fit for legged robots working on the Earth, such as BigDog [11] and HyQ2Max [12], an electric actuator is considered to be more realistic to use on the Moon. For example, the hexapod SpaceClimber [13], whose legs driven by 24 motors were arranged symmetrically along the sagittal plane of the body, was designed to execute extraterrestrial exploration. This probe succeeded to provide mobility in slopes with inclinations of up to 40° and capacity of negotiating obstacles of up to 40 cm height. An eight-legged scorpion driven by 32 motors [14] was further designed by this group to enhance the terrain adaptability. Electric quadruped robots, such as MIT Cheetah 3 [15], ANYmal [16], and SpaceBok [17,18], are more widely studied, owing to their less driving modules. The SpaceBok was designed for roving on the Moon; it used parallel elastic elements to store and release energy and mounted a reaction wheel to adjust the pitch angle for achieving continuous jumping gait. These legged robots show good agility and traversability in complex terrain. However, the collision between body and ground is prone to occur during landing due to the short buffer stroke and low deceleration ratio in leg design, which will damage the recoil engine nozzle installed at the bottom of the robot.

Movable landers, which can perform a buffer landing and complete walking on the lunar surface, are developed [19]. A legged mobile lander using four identical 2-UPS&UP parallel mechanism was designed by Lin et al. [20], which integrated the functions of landing, walking, deploying, and orientation adjusting. To improve the landing impact performance of a movable lander, Zhou et al. [21] proposed a newly designed gradual energy-absorbing structure subjected to impact loads. Lin et al. [22,23] addressed a novel synthesis method for legged mobile landers using truss-mechanism transformation based on the structure of Chang’e lander. Han et al. [24,25] reported the dimensional design and synthesis methods for cantilever-type legged lander using performance atlas. All these mobile landers have not completed the prototype design, and their buffer landing performance needs to be further tested. The walking function of these mobile landers is achieved by adding an additional motor to each leg branch chain, with impact absorbed by internal aluminum honeycomb material during landing, which leaves the problem of irreversible deformation of the buffer material unsolved. Thus, these landers are still only capable of one-off buffer landing.

The above two traditional strategies to improve the probe performance still have some limitations. Scientists have proposed a novel method to achieve long-range mobility by flying on a planet with a thin atmosphere. The conceptual design of Dragonfly lander [26] was proposed to achieve landing on the Titan, in which four rotors were used and many instruments, including sampling drill, cameras, and multimission radioisotope thermoelectric, were carried. This conceptual design was proved to be feasible as the Ingenuity helicopter [27], with a weight of 1.8 kg, a height of 0.5 m, and a peak motor power of 350 W, successfully achieved multiple takeoffs and landings on Mars. This flying approach can be borrowed to the Moon exploration. As shown in Fig.1, the lander utilizes the recoil engines to hover in the air and then uses the vision system to find a relatively flat terrain for buffer landing. After completing the exploration task near the current landing site as a rover, the recoil engines will rework and drive the lander to fly over complex terrains, such as valleys or mountains, and the above procedures will be repeated to achieve buffer landing on another site. On the basis of this assumption, a repetitive hexapod robot integrated with walking and landing functions is designed by our team. However, this hexapod robot has more driving modules and larger real-time communication data with controller compared with the four-legged lander. These features require more powerful recoil engines and more fuel, resulting in a lethal problem of higher error rate, so the detailed landing ability in the cases of leg failure was systematically discussed in our previous study [28].

Fig.1 Conceptual sketch of repetitive landing.

Full size|PPT slide

In this study, a four-legged landing and walking robot called “FLLWR” is designed for lunar exploration. It integrates the functions of repetitive lander and rover on a lightweight body and has less communication data simultaneously. The impact energy is converted into heat energy on motors through hybrid compliance to avoid the damage of internal structure, in which a passive spring is mounted between the thigh and shank to provide passive compliance, and three integrated drive units (IDUs) are mounted on each leg to provide active compliance. During the landing period, each IDU imitates the dynamics of an active torsional spring and an active torsional damper by a compliance control method [29,30]. Our robot is promising to achieve repetitive flying and landing by recoil engines because this approach has protected the intact functions of legs, especially in challenging terrains, such as deep craters or mountaintops on the Moon [31]. The landing control method is a critical technical challenge for the FLLWR because it determines the success of landing. A novel method with landing method framework, state switching strategy, trajectory planning, compliance planning, and height and attitude stabilizer, is proposed to solve this problem. The results of numerical simulations and experiments on a five-degree-of-freedom lunar gravity testing platform (5-DoF-LGTP) without/with horizontal disturbance velocity simultaneously show that this novel method can be used to achieve buffer landing for the FLLWR.

The rest of this paper is organized as follows. Section 2 introduces the robot design. Section 3 establishes the kinematics and dynamics of quadruped robots. Section 4 proposes the landing control method. Section 5 gives the landing numerical simulation. Section 6 discusses the experimental results on the 5-DoF-LGTP. Section 7 provides the conclusions.

2 2 Four-legged landing and walking robot

As shown in Fig.2, an FLLWR is designed to implement three functions: actively controlled landing, roving exploration, and repetitive landing. Its size is about 1.35 m length, 1.35 m width, and 0.70 m height, and its weight is about 40 kg. The robot mainly consists of one body and four identical legs. All legs are evenly and symmetrically distributed under the body. The body mainly carries the control system, sensing system, power systems, and other payloads. The control system consists of motor drivers, I/O modules, and motion controller. The sensing system includes inertial measurement unit (IMU), laser radar, encoders, torque sensors, and others. The power system includes battery packages. Load weights are used to emulate other payloads for lunar exploration. The body is covered by a square pyramid shell made of aluminum alloy, and the internal electronics are protected by the shell.

Fig.2 Structure of a four-legged landing and walking robot.

Full size|PPT slide

2.1 2.1 Multifunction design

The robot has multiple functions, including folding, unfolding, repetitive landing, and walking. The demonstrations for these functions are given in Fig.3. The robot keeps folding morphology during the launching period on the Earth with all legs folded under the body, and this compact structure reduces the occupied volume in rocket. When the robot reaches the predetermined landing site, the propulsion engine begins to decelerate and stabilize the attitude of the robot until the landing conditions are met. When the robot is about to land on the Moon, the unfolding process is initiated to extend all legs for lengthening compression stroke and enhancing support stability. The propulsion engine is turned off, the robot begins free falling, and the landing control method starts work at the same time. Each leg is controlled in accordance with the foot-ground touch events and the foot-ground contact forces.

Fig.3 Multifunction design during different phase: (a) buffer landing like a cantilever lander and (b) roving like a four-legged robot.

Full size|PPT slide

The robot becomes a four-legged rover adopting a stable wave gait [32] with at most one leg swinging and at least three legs supporting after successful landing. The walking process for one gait period is shown in Fig.3(b). In locomotion, the green arrow line indicates the center of mass (COM) projection movement from the old green point to the new red point. The red arrow arc indicates that one leg swings, and the triangle formed by blue dotted line indicates the supporting polygon.

2.2 2.2 Leg design

The prototype and mechanism of the leg are presented in Fig.4. The leg mechanism is a hybrid parallel-serial mechanism, where a parallel five-bar mechanism (PFBM) is serially connected with a revolute pair (RP). The PFBM is a parallelogram planar parallel mechanism, where the connecting link is parallel with the thigh, and the rocker arm is parallel with the shank. Each leg has three active degrees-of-freedom (DoFs) and is actuated by three IDUs. The PFBM is used to realize the swing and stretching/shrinking movement of the leg in the leg sagittal plane. The RP between the body and the PFBM is used to realize the abduction/adduction movement of the leg in the body longitudinal plane. A passive linear spring is added between the thigh and the shank to balance the partial weight of the lander and reduce the required driving torque and the impact on robot from foot-ground event.

Fig.4 Three active degrees-of-freedom hybrid leg design: (a) leg prototype and (b) leg mechanism. IDU: integrated drive unit, RP: revolute pair, PFBM: parallel five-bar mechanism.

Full size|PPT slide

2.3 2.3 Actuator design

An external rotor motor rather than an internal one is selected to provide more powerful drive capability in the same volume. A reducer is required in the actuator to provide the high torque during landing. The motor and the reducer are then integrated and designed as an IDU, consisting of servo motor, reducer, torque sensor, encoder, brake, bearing, coupler, and shell. Fig.5 shows the prototype of IDU and the main components. In accordance with the torque and velocity requirements, the motor parameters are selected as follows: peak torque of 2 N·m and maximum speed of 5000 r/min, the reducer parameters are selected as transmission ratio of 100 and maximum output torque of 200 N·m. For the feedback control, the rotary encoder with 4096 counts per turn is selected. The torque sensor has a precision of 0.1% and a measurement range of ±250 N·m. The IDU is compact in structure and weighs 1.1 kg.

Fig.5 Design of the integrated drive unit: (a) three-dimensional and (b) inner structure.

Full size|PPT slide

2.4 2.4 Control system

The hardware structure of the control system is illustrated in Fig.6. The control algorithm runs on a real-time operating system in the controller. The roll and pitch angles of the body are sensed by IMU at 0.5 kHz. On the basis of the EtherCAT bus system at 4 kHz, the actual torque ${\mathit{\tau }}_{\mathrm{a}\mathrm{j}}$, actual position ${\mathit{\theta }}_{\mathrm{a}\mathrm{j}}$, and actual velocity ${\dot{\mathit{\theta }}}_{\mathrm{a}\mathrm{j}}$ of each joint are detected by the torque sensor and encoder. Remote communication is used to achieve the human–robot interaction. The 12 command joint torques ${\mathit{\tau }}_{\mathrm{c}\mathrm{j}}$ are transmitted to the IDUs by the EtherCAT bus as well.

Fig.6 Hardware structure of motion controller. RTOS: real time operating system.

Full size|PPT slide

3 3 Kinematics and dynamics of quadruped robots

3.1 3.1 Coordinate frames

Several coordinate frames are defined to establish the robot model, as shown in Fig.7. For the body coordinate frame ${\mathrm{\Sigma }}_{\mathrm{b}}$, the origin $O_{\mathrm{b}}$ is located at the body center, $x_{\mathrm{b}}$ is the forward direction, $y_{\mathrm{b}}$ is the left direction, and ${\mathit{\text{z}}}_{\mathrm{b}}$ is the upper direction. For the world coordinate frame ${\mathrm{\Sigma }}_{\mathrm{w}}$, its directions coincide with ${\mathrm{\Sigma }}_{\mathrm{b}}$ at the initial moment, and the origin $O_{\mathrm{w}}$ is located at the projection point of $O_{\mathrm{b}}$ on the ground. For the ith leg coordinate frame ${\mathrm{\Sigma }}_{\mathrm{l}i}$, the origin $O_{\mathrm{l}i}$ is located at the center point of the hip joint, $x_{\mathrm{l}i}$ points to the body center, and $y_{\mathrm{l}i}$ is parallel to ${\mathit{\text{z}}}_{\mathrm{b}}$. For the IMU coordinate frame ${\mathrm{\Sigma }}_{\mathrm{i}\mathrm{m}\mathrm{u}}$, the axis directions are the same as those of ${\mathrm{\Sigma }}_{\mathrm{b}}$.

Fig.7 Coordinate system in motion planning: (a) robot coordinate and (b) leg coordinate.

Full size|PPT slide

3.2 3.2 Leg modeling

The leg model includes leg forward kinematic model, leg inverse kinematic model, velocity model, and force model. The leg forward kinematic model is used to establish the relationship between the joint angles and the tiptoe position in the leg coordinate frame. The joint angles form a generalized coordinate vector $\mathit{q}=[\begin{array}{ccc}{\theta }_{\mathrm{s}}& {\theta }_{\mathrm{t}}& {\theta }_{\mathrm{r}}\end{array}]{}^{\mathrm{T}}$, where ${\theta }_{\mathrm{s}}$ is the side angle, ${\theta }_{\mathrm{t}}$ is the thigh angle, and ${\theta }_{\mathrm{r}}$ is the rocker arm angle. The tiptoe position ${\mathit{P}}_{\mathrm{t}\mathrm{i}\mathrm{p}}=[\begin{array}{ccc}{x}_{\mathrm{t}\mathrm{i}\mathrm{p}}& y_{\mathrm{t}\mathrm{i}\mathrm{p}}& {\mathit{\text{z}}}_{\mathrm{t}\mathrm{i}\mathrm{p}}\end{array}{]}^{\mathrm{T}}$ can be obtained as follows:

(1)${\mathit{P}}_{\mathrm{t}\mathrm{i}\mathrm{p}}=\left[\begin{array}{c}{L}_{\mathrm{t}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{t}}-L_{\mathrm{s}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{r}}\\ \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{s}}\left(L_{\mathrm{t}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{t}}-L_{\mathrm{s}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{r}}\right)\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{s}}\left(L_{\mathrm{t}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{t}}-L_{\mathrm{s}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{r}}\right)\end{array}\right],$

where $L_{\mathrm{t}}$ and $L_{\mathrm{s}}$ are the lengths of thigh and shank, respectively.

The leg inverse kinematic model is employed to calculate the joint angles from the tiptoe position in the leg coordinate frame. For the given tiptoe position, we have

(2)${\theta }_{\mathrm{s}}=\arctan \frac{{\mathit{\text{z}}}_{\mathrm{t}\mathrm{i}\mathrm{p}}}{y_{\mathrm{t}\mathrm{i}\mathrm{p}}},$

(3)${\theta }_{\mathrm{t}}=\arcsin \frac{-L_{\mathrm{s}}\sin B}{\sqrt{x_{\mathrm{t}\mathrm{i}\mathrm{p}}^2+y_{\mathrm{t}\mathrm{i}\mathrm{p}}^2/\phantom{y_{\mathrm{t}\mathrm{i}\mathrm{p}}^2{\cos }^2{\theta }_{\mathrm{s}}}{\cos }^2{\theta }_{\mathrm{s}}}}+\arctan \frac{y_{\mathrm{t}\mathrm{i}\mathrm{p}}}{x_{\mathrm{t}\mathrm{i}\mathrm{p}}\cos {\theta }_{\mathrm{s}}},$

(4)${\theta }_{\mathrm{r}}={\theta }_{\mathrm{t}}-B,$

where

(5)$B=\arccos \frac{-x_{\mathrm{t}\mathrm{i}\mathrm{p}}^2-y_{\mathrm{t}\mathrm{i}\mathrm{p}}^2/\phantom{-x_{\mathrm{t}\mathrm{i}\mathrm{p}}^2-y_{\mathrm{t}\mathrm{i}\mathrm{p}}^2{\cos }^2{\theta }_{\mathrm{s}}+L_{\mathrm{t}}^2+L_{\mathrm{s}}^2}{\cos }^2{\theta }_{\mathrm{s}}+L_{\mathrm{t}}^2+L_{\mathrm{s}}^2}{2L_{\mathrm{t}}{L}_{\mathrm{s}}}.$

In accordance with Eqs. (2)‒(4), the generalized vector of the input angles $\mathit{q}=[\begin{array}{ccc}{\theta }_{\mathrm{s}}& {\theta }_{\mathrm{t}}& {\theta }_{\mathrm{r}}\end{array}]{}^{\mathrm{T}}$ is obtained.

Let $\dot{\mathit{q}}=[\begin{array}{ccc}{\dot{\theta }}_{\mathrm{s}}& {\dot{\theta }}_{\mathrm{t}}& {\dot{\theta }}_{\mathrm{r}}\end{array}{]}^{\mathrm{T}}$ be the generalized vector of the joint velocity. The tiptoe velocity ${\dot{\mathit{P}}}_{\mathrm{t}\mathrm{i}\mathrm{p}}=[\begin{array}{ccc}{\dot{x}}_{\mathrm{t}\mathrm{i}\mathrm{p}}& {\dot{y}}_{\mathrm{t}\mathrm{i}\mathrm{p}}& {\dot{\mathit{\text{z}}}}_{\mathrm{t}\mathrm{i}\mathrm{p}}\end{array}]{}^{\mathrm{T}}$ by differentiating Eq. (1) can be expressed as

(6)${\dot{\mathit{P}}}_{\mathrm{t}\mathrm{i}\mathrm{p}}={\mathit{J}}_{\mathrm{v}}(\mathit{q})\dot{\mathit{q}},$

where ${\mathit{J}}_{\mathrm{v}}(\mathit{q})$ is the velocity Jacobian matrix as follows:

(7)$\begin{array}{rl}& {\mathit{J}}_{\mathrm{v}}(\mathit{q})=\\ & \left[\begin{array}{ccc}0& -L_{\mathrm{t}}\sin {\theta }_{\mathrm{t}}& L_{\mathrm{s}}\sin {\theta }_{\mathrm{r}}\\ -\sin {\theta }_{\mathrm{s}}\left(L_{\mathrm{t}}\sin {\theta }_{\mathrm{t}}-L_{\mathrm{s}}\sin {\theta }_{\mathrm{r}}\right)& L_{\mathrm{t}}\cos {\theta }_{\mathrm{s}}\cos {\theta }_{\mathrm{t}}& -L_{\mathrm{s}}\cos {\theta }_{\mathrm{s}}\cos {\theta }_{\mathrm{r}}\\ \cos {\theta }_{\mathrm{s}}\left(L_{\mathrm{t}}\sin {\theta }_{\mathrm{t}}-L_{\mathrm{s}}\sin {\theta }_{\mathrm{r}}\right)& L_{\mathrm{t}}\sin {\theta }_{\mathrm{s}}\cos {\theta }_{\mathrm{t}}& -L_{\mathrm{s}}\sin {\theta }_{\mathrm{s}}\cos {\theta }_{\mathrm{r}}\end{array}\right].\end{array}$

The length of passive spring can be written as follows:

(8)$L_{\mathrm{p}\mathrm{s}}=\sqrt{L_{\mathrm{r}}^2+L_{\mathrm{f}}^2-2L_{\mathrm{l}}{L}_{\mathrm{r}}\cos (\text{π}-{\theta }_{\mathrm{t}}+{\theta }_{\mathrm{r}})},$

where $L_{\mathrm{l}}$ and $L_{\mathrm{r}}$ are the distances from the intersection point of thigh and shank to the left and right fixed points of spring, respectively, and $L_{\mathrm{p}\mathrm{s}}$ is the current length of the passive spring.

Its elastic potential energy of passive spring ($E_{\mathrm{p}\mathrm{s}}$) is

(9)$E_{\mathrm{p}\mathrm{s}}=\frac{K_{\mathrm{p}\mathrm{s}}{\left(L_{\mathrm{p}\mathrm{s}}-L_{\mathrm{p}\mathrm{s}0}\right)}^2}{2},$

where $K_{\mathrm{p}\mathrm{s}}$ and $L_{\mathrm{p}\mathrm{s}0}$ are the stiffness coefficient and original length of the passive spring, respectively.

The virtual work of passive spring ($\text{δ}{W}_{\mathrm{p}\mathrm{s}}$) is

(10)$\text{δ}{W}_{\mathrm{p}\mathrm{s}}=-\frac{\mathrm{\partial }{E}_{\mathrm{p}\mathrm{s}}}{\mathrm{\partial }\mathit{q}}\text{δ}\mathit{q},$

where $\text{δ}\mathit{q}$ is the virtual displacement of the generalized coordinate vector $\mathit{q}$.

The relationship between the tiptoe force ${\mathit{F}}_{\mathrm{t}\mathrm{i}\mathrm{p}}=$ $[\begin{array}{ccc}{F}_x& F_y& F_{\mathit{\text{z}}}\end{array}{]}^{\mathrm{T}}$ and the joint torque vector $\mathit{\tau }=[\begin{array}{ccc}{\tau }_{\mathrm{s}}& {\tau }_{\mathrm{t}}& {\tau }_{\mathrm{r}}\end{array}{]}^{\mathrm{T}}$ can be obtained by using the d’Alembert virtual work principle as follows:

(11)${\mathit{F}}_{\mathrm{t}\mathrm{i}\mathrm{p}}^{\mathrm{T}}\text{δ}{\mathit{P}}_{\mathrm{t}\mathrm{i}\mathrm{p}}=\left({\mathit{\tau }}^{\text{T}}-{\mathit{H}}^{\text{T}}\right)\text{δ}\mathit{q},$

(12)${\mathit{H}}^{\mathrm{T}}=[\begin{array}{ccc}0& h& -h\end{array}],$

(13)$h=\frac{-K_{\mathrm{p}\mathrm{s}}\left(L_{\mathrm{p}\mathrm{s}}-L_{\mathrm{p}\mathrm{s}0}\right)L_{\mathrm{l}}{L}_{\mathrm{r}}\sin ({\theta }_{\mathrm{t}}-{\theta }_{\mathrm{r}})}{L_{\mathrm{p}\mathrm{s}}},$

where $\text{δ}{\mathit{P}}_{\mathrm{t}\mathrm{i}\mathrm{p}}$ is the virtual displacements of tiptoe and active joint, and H is the correction vector.

Considering the derivative relationship in Eq. (6), Eq. (11) can be rewritten as follows:

(14)$\mathit{\tau }={\mathit{J}}_{\mathrm{v}}(\mathit{q}{)}^{\mathrm{T}}{\mathit{F}}_{\mathrm{t}\mathrm{i}\mathrm{p}}+\mathit{H}.$

An additional correction vector H is added to the joint torque vector τ due to the use of passive spring. When the leg is compressed for buffering, the passive spring reduces the thigh torque amplitude and increases the shank torque amplitude. The peak torque is reduced, and the buffering capacity is improved because the amplitude of thigh torque is greater than that of shank torque in the buffering process.

3.3 3.3 Body modeling

The body model includes kinematic model and dynamic model. The kinematic model is used to express the relationship of position vector between the leg coordinate frame ${\mathrm{\Sigma }}_{\mathrm{l}i}$ and the world coordinate frame ${\mathrm{\Sigma }}_{\mathrm{w}}$. As shown in Fig.7, ${}^{\mathrm{l}i}{\mathit{P}}_{\mathrm{t}\mathrm{i}\mathrm{p}}$ is the tiptoe position in the leg coordinate frame ${\mathrm{\Sigma }}_{\mathrm{l}i}$. ${}^{\mathrm{b}}{\mathit{P}}_{\mathrm{l}i}$ is used to denote the origin position $O_{\mathrm{l}i}$ of ${\mathrm{\Sigma }}_{\mathrm{l}i}$ in the body coordinate frame ${\mathrm{\Sigma }}_{\mathrm{b}}$. ${}^{\mathrm{b}}{\mathit{R}}_{\mathrm{l}i}$ is the rotation transformation matrix from ${\mathrm{\Sigma }}_{\mathrm{l}i}$ to ${\mathrm{\Sigma }}_{\mathrm{b}}$. ${}^{\mathrm{w}}{\mathit{P}}_{\mathrm{b}}$ is the body position in the world coordinate frame ${\mathrm{\Sigma }}_{\mathrm{w}}$. ${}^{\mathrm{w}}{\mathit{R}}_{\mathrm{b}}$ is the rotation transformation matrix from ${\mathrm{\Sigma }}_{\mathrm{b}}$ to ${\mathrm{\Sigma }}_{\mathrm{w}}$. The tiptoe position ${}^{\mathrm{w}}{\mathit{P}}_{\mathrm{t}\mathrm{i}\mathrm{p}}$ in the world coordinate frame is given as follows:

(15)${}^{\mathrm{w}}{\mathit{P}}_{\mathrm{t}\mathrm{i}\mathrm{p}}={}^{\mathrm{w}}{\mathit{P}}_{\mathrm{b}}+{}^{\mathrm{w}}{\mathit{R}}_{\mathrm{b}}\left({}^{\mathrm{b}}{\mathit{P}}_{\mathrm{l}i}+{}^{\mathrm{b}}{\mathit{R}}_{\mathrm{l}i}{}^{\mathrm{l}i}{\mathit{P}}_{\mathrm{t}\mathrm{i}\mathrm{p}}\right).$

A virtual dynamic model is established for the gravity compensation and the attitude adjustment. The mass of the leg can be ignored because it is smaller relative to the mass of the body to simplify the dynamics model. The contact between the tiptoe and the ground is inelastic, and no slipping or rebounding occurs. As shown in Fig.7, ${\mathit{r}}_{\mathrm{c}\mathrm{o}\mathrm{m}}$ is the vector from the COM to the origin $O_{\mathrm{b}}$. ${\mathit{r}}_i$ is the vector from the origin $O_{\mathrm{b}}$ to the ith tiptoe. ${\mathit{F}}_i$ is the ith tiptoe force from the ground. The body dynamics can be formulated as

(16)$m_{\mathrm{b}}{\mathit{a}}_{\mathrm{b}}=m_{\mathrm{b}}\mathit{g}+\sum _{i=1}^4{\mathit{F}}_i,$

(17)$\frac{\mathrm{d}}{\mathrm{d}t}\left({\mathit{I}}_{\mathrm{b}}{\text{ω}}_{\mathrm{b}}\right)=\sum _{i=1}^4\left(\left({\mathit{r}}_{\mathrm{c}\mathrm{o}\mathrm{m}}+{\mathit{r}}_i\right)\times {\mathit{F}}_i\right),$

where $t$, $m_{\mathrm{b}}$, ${\mathit{I}}_{\mathrm{b}}$, ${\mathit{a}}_{\mathrm{b}}$, ${\text{ω}}_{\mathrm{b}}$, and $\mathit{g}$ are the time, body mass, body inertia, body linear acceleration vector, body angular velocity vector, and gravitational acceleration vector, respectively.

4 4 Landing control method

4.1 4.1 Landing method framework

The landing method controls the FLLWR to absorb the impact energy in the leg compression stroke to realize buffer landing. The maximum torque of IDU can be controlled to be less than the peak torque. Four legs are coordinated in accordance with the states of the legs.

As shown in Fig.8, the landing method framework consists of three layers: hardware input layer, landing method layer, and hardware output layer. The hardware input layer processes the actual joint torque ${\mathit{\tau }}_{\mathrm{a}\mathrm{j}}$, the actual joint angles ${\mathit{\theta }}_{\mathrm{a}\mathrm{j}}$, the actual roll angle ${\theta }_{\mathrm{a}\mathrm{r}}$, the actual pitch angle ${\theta }_{\mathrm{a}\mathrm{p}}$, and the corresponding velocities ${\dot{\mathit{\theta }}}_{\mathrm{a}\mathrm{j}}$, ${\dot{\theta }}_{\mathrm{a}\mathrm{r}}$, and ${\dot{\theta }}_{\mathrm{a}\mathrm{p}}$ from the torque sensors, encoders, and IMU. In the landing method layer, if the amplitude of the tiptoe force for each leg in the vertical direction exceeds a threshold during the landing period, which signifies the occurrence of a touching-ground event, then the leg state will be changed to the grounded phase from the aerial phase. The trajectory planning module gives the reference trajectory of the joint angle ${\mathit{\theta }}_{\mathrm{r}\mathrm{j}}$, the reference roll angle ${\theta }_{\mathrm{r}\mathrm{r}}$, the reference pitch angle ${\theta }_{\mathrm{r}\mathrm{p}}$, and the corresponding velocities ${\dot{\mathit{\theta }}}_{\mathrm{r}\mathrm{j}}$, ${\dot{\theta }}_{\mathrm{r}\mathrm{r}}$, and ${\dot{\theta }}_{\mathrm{r}\mathrm{p}}$. Simultaneously, the active stiffness coefficient K and the damping coefficient B are given by the compliance planning. The actual tiptoe position can be obtained in real time by using the forward kinematics. The stabilizer is used to control the body roll angle, the body pitch angle, and the body height by generating the reference joint torque ${\mathit{\tau }}_{\mathrm{r}\mathrm{j}}$. The 12 joint impedance controllers generate the command joint torque ${\mathit{\tau }}_{\mathrm{c}\mathrm{j}}$ of each IDU by using Eq. (18):

Fig.8 Landing method framework consisting of three layers: hardware inputs, landing method, and hardware outputs.

Full size|PPT slide

(18)${\mathit{\tau }}_{\mathrm{c}\mathrm{j}}-{\mathit{\tau }}_{\mathrm{r}\mathrm{j}}=K\left({\mathit{\theta }}_{\mathrm{r}\mathrm{j}}-{\mathit{\theta }}_{\mathrm{a}\mathrm{j}}\right)+B\left({\dot{\mathit{\theta }}}_{\mathrm{r}\mathrm{j}}-{\dot{\mathit{\theta }}}_{\mathrm{a}\mathrm{j}}\right).$

In the hardware output layer, the drivers transform ${\mathit{\tau }}_{\mathrm{c}\mathrm{j}}$ into the corresponding current to control the motors.

4.2 4.2 State switching strategy

As illustrated in Fig.9, the robot state consists of the leg state and the body state. The leg state includes aerial phase and grounded phase, and the body state includes landing preparation phase, synergy adjustment phase, and full-supporting phase. The robot reaches above the landing site and uses the engines to adjusts the velocity and attitude to meet the landing condition. The lander starts to descend under the lunar gravity after turning off the engines. Concurrently, each leg switches to the aerial phase, and the body switches to the landing preparation phase. The legs touch the ground successively owing to the irregular terrain after a period of descent. Making a quick response to touch-ground event by the actual tiptoe force is difficult because its change is relatively small owing to the low coefficients of stiffness and damping in the landing prepare phase (Section 4.4). The virtual tiptoe force ${\mathit{F}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$ generated by high virtual coefficients of stiffness ($K_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$) and damping ($B_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$) replaces the actual tiptoe force to judge the touchdown and to increase the response sensitivity to touching-ground event. ${\mathit{F}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$ can be calculated as follows:

Fig.9 State switching method during buffer landing.

Full size|PPT slide

(19)${\mathit{\tau }}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}=K_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}\left({\mathit{\theta }}_{\mathrm{r}\mathrm{j}}-{\mathit{\theta }}_{\mathrm{a}\mathrm{j}}\right)+B_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}\left({\dot{\mathit{\theta }}}_{\mathrm{r}\mathrm{j}}-{\dot{\mathit{\theta }}}_{\mathrm{a}\mathrm{j}}\right),$

(20)${\mathit{F}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}={\left({\mathit{J}}_{\mathrm{v}}{(\mathit{q})}^{\mathrm{T}}\right)}^{-1}\left({\mathit{\tau }}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}-\mathit{H}\right),$

where ${\mathit{\tau }}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$ is the virtual joint torque.

When the touching-ground event occurs at the first time, the body switches to the synergy adjustment phase in which some legs are grounded phase and the others are aerial phase. When all legs switch to the grounded phase, the robot state switches to the full-supporting phase.

4.3 4.3 Trajectory planning

The trajectory planning scheme for leg 4 is illustrated in Fig.10, and the other legs also adopt the same planning. In the landing preparation phase, the world coordinate system is assumed to coincide exactly with the body coordinate system. Therefore, the body position in the world coordinate system is zero, and the tiptoes keep initial positions: $x_{\mathrm{t}\mathrm{i}\mathrm{p}0}$, $y_{\mathrm{t}\mathrm{i}\mathrm{p}0}$, and ${\mathit{\text{z}}}_{\mathrm{t}\mathrm{i}\mathrm{p}0}$. An active compression distance of $d$ is planned to reduce the impact because the impact energy is mainly distributed in the vertical direction. The initial tiptoe positions of all legs used in the experiments are listed in Tab.1. In the synergy adjustment phase, if the leg state changes to the grounded phase from the aerial phase, the net output torques of the thigh and shank IDUs will become 0 N·m due to the zero tiptoe force state (Section 4.4), resulting in a passive shrink under the action of the downward body. This design has two advantages. The lengths of the grounded legs will become shorter continuously, and the lengths of the aerial legs retain their initial value, which will contribute to the good adaption to landing terrain. The grounded legs will exert minimal influence on the body position and attitude because their forces are close to zero. At the beginning of the full-supporting phase, the world coordinate frame is updated to coincide with the current body coordinate frame, and the actual tiptoe positions, including $x_{\mathrm{t}\mathrm{i}\mathrm{p}1}$, $y_{\mathrm{t}\mathrm{i}\mathrm{p}1}$, and ${\mathit{\text{z}}}_{\mathrm{t}\mathrm{i}\mathrm{p}1}$, are recorded as a referenced tiptoe position in the whole supporting phases. For the body, an active body compression in the z direction from initial position ${\mathit{\text{z}}}_{\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}0}$ to target position ${\mathit{\text{z}}}_{\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}1}$ is planned to reduce the peak torque. The real-time body interpolation trajectory ${\mathit{\text{z}}}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{p}}$ can be obtained from the following equation:

Tab.1 Initial tiptoe positions before buffer landing

Leg	xtip0/m	ytip0/m	ztip0/m
1	0.788	0.788	−0.750
2	−0.788	0.788	−0.750
3	−0.788	−0.788	−0.750
4	0.788	−0.788	−0.750

Fig.10 Trajectory planning of the body and tiptoe of leg 4 during buffer landing.

Full size|PPT slide

(21)${\mathit{\text{z}}}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{p}}=c_0+c_1\left(t-t_0\right)+c_2{\left(t-t_0\right)}^2+c_3{\left(t-t_0\right)}^3+c_4{\left(t-t_0\right)}^4+c_5{\left(t-t_0\right)}^5+c_6{\left(t-t_0\right)}^6+c_7{\left(t-t_0\right)}^7,$

(22)$\{\begin{array}{l}[\begin{array}{cccc}{c}_0& c_1& c_2& c_3\end{array}]=[\begin{array}{cccc}{\mathit{\text{z}}}_{\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}0}& v_0& a_0/2& j_0/6\end{array}],\\ c_4=\left(210d-T\left(\left(30a_0-15a_1\right)T+\left(4j_0+j_1\right)T^2+120v_0+90v_1\right)\right)/\left(6T^4\right),\\ c_5=\left(-168d+T\left(\left(20a_0-14a_1\right)T+\left(2j_0+j_1\right)T^2+90v_0+78v_1\right)\right)/\left(2T^5\right),\\ c_6=\left(420d-T\left(\left(20a_0-39a_1\right)T+\left(4j_0+3j_1\right)T^2+216v_0+204v_1\right)\right)/\left(6T^6\right),\\ c_7=\left(-120d+T\left(\left(12a_0-12a_1\right)T+\left(j_0+j_1\right)T^2+60v_0+60v_1\right)\right)/\left(6T^7\right),\\ [\begin{array}{cc}T& d\end{array}]=[\begin{array}{cc}{t}_1-t_0& {\mathit{\text{z}}}_{\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}1}-{\mathit{\text{z}}}_{\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}0}\end{array}],\end{array}$

where $v_0$, $v_1$, $a_0$, $a_1$, $j_0$, and $j_1$ are the initial velocity, target velocity, initial acceleration, target acceleration, initial jerk, and target jerk of body in the z direction, respectively, $t_0$, $t_1$, $t$, and $T$ are the initial time, end time, current time, and total duration time, respectively, and $c_i$ $(i=1,2,...,7)$ is the coefficient of interpolation trajectory. The equation can produce a trajectory of flexible impact by setting the velocity, acceleration, and jerk at the beginning and end time to zero. It avoids computing the inverse of 8 by 8 matrix in real time. The angle of each joint can be obtained by the inverse kinematics after completing the trajectory planning of the body and tiptoes during the whole landing period.

4.4 4.4 Compliance planning

The mechanical property of the leg is determined by the stiffness and damping coefficients, which can be adjusted actively through the compliance control. The compliance planning method for thigh and shank joints is shown in Fig.11, and the side joint keeps the high coefficients of stiffness ($K_2$) and damping ($B_3$) all the time. All thigh and shank joints keep the low stiffness coefficients $K_0$ and low damping coefficients $B_0$ to avoid a huge impact on each leg at the touching-ground moment in the landing preparation phase. At the beginning of synergy adjustment phase, the thigh and shank joints of grounded legs smoothly switch to zero tiptoe force state by the transition 1, and those of the aerial legs still retain the initial values $K_0$ and $B_0$. The damping coefficients increase to a median value $B_2$ by the transition 2 to avoid the torque overload because the angular velocity of active joint may reach its maximum at the start of full-supporting phase. Given that the angular velocity decreases after the energy dissipation in the transition 2, the damping coefficient can reach the maximum $B_3$ by the transition 3 to release the motor torque potential and realize the energy quick dissipation. The stiffness coefficient reaches the maximum $K_2$ by the transition 2 because the elastic force increases slowly from zero without an amplitude step. The switch equation of stiffness and damping coefficients is written as follows:

Fig.11 Compliance planning of stiffness and damping coefficients for thigh and shank joints.

Full size|PPT slide

(23)$\begin{array}{rl}co{e}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{p}}=& -3c_0{\left(1-t_{\mathrm{r}}\right)}^2-6c_1\left(1-t_{\mathrm{r}}\right)t_{\mathrm{r}}+3c_1{\left(1-t_{\mathrm{r}}\right)}^2\\ & -3c_2{t}_{\mathrm{r}}^2+6c_2\left(1-t_{\mathrm{r}}\right)t_{\mathrm{r}}+3c_3{t}_{\mathrm{r}}^2,\end{array}$

(24)$\{\begin{array}{l}[\begin{array}{cc}{c}_0& c_3\end{array}]=[\begin{array}{cc}co{e}_0& co{e}_1\end{array}],\\ c_1=Tcoe{V}_0/3+co{e}_0,\\ c_2=-Tcoe{V}_1/3+co{e}_1,\\ T=t_1-t_0,\\ t_{\mathrm{r}}=\left(t-t_0\right)/T,\end{array}$

where $co{e}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{p}}$, $co{e}_0$, $co{e}_1$, $coe{V}_0$, $coe{V}_1$, and $t_{\mathrm{r}}$ are the current interpolation coefficient, initial coefficient, target coefficient, initial coefficient velocity, target coefficient velocity, and duration ratio time (from 0 to 1), respectively. The equation avoids solving the inverse of 4 by 4 matrix due to the explicit expression of $c_0$, $c_1$, $c_2$, and $c_3$.

4.5 4.5 Height and attitude stabilizer design

The resultant force in the vertical direction, including the leg force, the gravity force, and the inertial force, should be kept zero all the time in the full-supporting phase to maintain a stable height. When the roll and pitch angles of the body are not equal to the reference value, an extra virtual torsion is needed to control the body attitude. In accordance with the dynamic model, we have

(25)$m_{\mathrm{b}}\mathit{g}+\sum _{i=1}^4{\mathit{F}}_{i\mathit{\text{z}}}=m_{\mathrm{b}}{\ddot{\mathit{z}}}_{\mathrm{b}},$

(26)$\sum _{i=1}^4\left(\left({\mathit{r}}_{\mathrm{c}\mathrm{o}\mathrm{m}}+{\mathit{r}}_i\right)\times {\mathit{F}}_{i\mathit{\text{z}}}\right)={\mathit{I}}_{\mathrm{b}}{\dot{\text{ω}}}_{\mathrm{b}},$

where ${\ddot{\mathit{z}}}_{\mathrm{b}}$ and ${\dot{\text{ω}}}_{\mathrm{b}}$ are the body linear and angular acceleration, and ${\mathit{F}}_{i\mathit{\text{z}}}$ is the ith leg force in vertical direction.

Here, the classical proportional–integral–derivative controllers are utilized to obtain the desired acceleration magnitude (${\ddot{\mathit{\text{z}}}}_{\mathrm{b}}$) and angular acceleration ${\dot{\text{ω}}}_{\mathrm{b}}=[\begin{array}{cc}{\dot{\omega }}_{\mathrm{b}x}& {\dot{\omega }}_{\mathrm{b}y}\end{array}{]}^{\mathrm{T}}$ as follows:

(27)${\ddot{\mathit{\text{z}}}}_{\mathrm{b}}=k_{\mathrm{p},\mathit{\text{z}}}\left({\mathit{\text{z}}}_{\mathrm{r}\mathrm{b}}-{\mathit{\text{z}}}_{\mathrm{a}\mathrm{b}}\right)+k_{\mathrm{i},\mathit{\text{z}}}\int \left({\mathit{\text{z}}}_{\mathrm{r}\mathrm{b}}-{\mathit{\text{z}}}_{\mathrm{a}\mathrm{b}}\right)\mathrm{d}t+k_{\mathrm{d},\mathit{\text{z}}}\left({\dot{\mathit{\text{z}}}}_{\mathrm{r}\mathrm{b}}-{\dot{\mathit{\text{z}}}}_{\mathrm{a}\mathrm{b}}\right),$

(28)$\left[\begin{array}{c}{\dot{\omega }}_{\mathrm{b}x}\\ {\dot{\omega }}_{\mathrm{b}y}\end{array}\right]=k_{\mathrm{p},\theta }\left[\begin{array}{c}{\theta }_{\mathrm{r}\mathrm{r}}-{\theta }_{\mathrm{a}\mathrm{r}}\\ {\theta }_{\mathrm{r}\mathrm{p}}-{\theta }_{\mathrm{a}\mathrm{p}}\end{array}\right]+k_{\mathrm{i},\theta }\left[\begin{array}{c}\int \left({\theta }_{\mathrm{r}\mathrm{r}}-{\theta }_{\mathrm{a}\mathrm{r}}\right)\mathrm{d}t\\ \int \left({\theta }_{\mathrm{r}\mathrm{p}}-{\theta }_{\mathrm{a}\mathrm{p}}\right)\mathrm{d}t\end{array}\right]+k_{\mathrm{d},\theta }\left[\begin{array}{c}{\dot{\theta }}_{\mathrm{r}\mathrm{r}}-{\dot{\theta }}_{\mathrm{a}\mathrm{r}}\\ {\dot{\theta }}_{\mathrm{r}\mathrm{p}}-{\dot{\theta }}_{\mathrm{a}\mathrm{p}}\end{array}\right],$

where $k_{\mathrm{p},\mathit{\text{z}}}$, $k_{\mathrm{i},\mathit{\text{z}}}$, and $k_{\mathrm{d},\mathit{\text{z}}}$ are the proportional, integral, and derivative gain of the body in vertical direction, respectively, ${\mathit{\text{z}}}_{\mathrm{r}\mathrm{b}}$, ${\mathit{\text{z}}}_{\mathrm{a}\mathrm{b}}$, ${\dot{\mathit{\text{z}}}}_{\mathrm{r}\mathrm{b}}$, and ${\dot{\mathit{\text{z}}}}_{\mathrm{a}\mathrm{b}}$ are the reference position, actual position, reference velocity, and actual velocity of the body in vertical direction, respectively, and $k_{\mathrm{p},\theta }$, $k_{\mathrm{i},\theta }$, and $k_{\mathrm{d},\theta }$ are the proportional, integral, and derivative gains of roll and pitch angles.

Rewriting Eqs. (25)‒(28) in matrix form yields

(29)$\begin{array}{rl}& \left[\begin{array}{cccc}1& 1& 1& 1\\ r_{\mathrm{c}\mathrm{o}\mathrm{m}y}+r_{1y}& r_{\mathrm{c}\mathrm{o}\mathrm{m}y}+r_{2y}& r_{\mathrm{c}\mathrm{o}\mathrm{m}y}+r_{3y}& r_{\mathrm{c}\mathrm{o}\mathrm{m}y}+r_{4y}\\ r_{\mathrm{c}\mathrm{o}\mathrm{m}x}+r_{1x}& r_{\mathrm{c}\mathrm{o}\mathrm{m}x}+r_{2x}& r_{\mathrm{c}\mathrm{o}\mathrm{m}x}+r_{3x}& r_{\mathrm{c}\mathrm{o}\mathrm{m}x}+r_{4x}\end{array}\right]\left[\begin{array}{c}{F}_{1\mathit{\text{z}}}\\ F_{2\mathit{\text{z}}}\\ F_{3\mathit{\text{z}}}\\ F_{4\mathit{\text{z}}}\end{array}\right]\\ =& \left[\begin{array}{c}{m}_{\mathrm{b}}\left({\ddot{\mathit{\text{z}}}}_{\mathrm{b}}-g\right)\\ I_{\mathrm{b}x}{\dot{\omega }}_{\mathrm{b}x}\\ I_{\mathrm{b}y}{\dot{\omega }}_{\mathrm{b}y}\end{array}\right].\end{array}$

where $g$ and $F_{iz}$ are the third component of $\mathit{g}$ and ${\mathit{F}}_{iz}$, respectively, $r_{\mathrm{c}\mathrm{o}\mathrm{m}x}$, $r_{\mathrm{c}\mathrm{o}\mathrm{m}y}$, $r_{ix}$, $r_{iy}$, $r_{iz}$, $I_{\mathrm{b}x}$, and $I_{\mathrm{b}y}$ are the corresponding components of ${\mathit{r}}_{\mathrm{c}\mathrm{o}\mathrm{m}}$, ${\mathit{r}}_i$, and ${\mathit{I}}_{\mathrm{b}}$.

Equation (29) is an indefinite equation and has many feasible solutions. Researchers have proposed many effective methods to compute the Moore–Penrose inverse matrix, and some other researchers have obtained the tiptoe forces by solving a quadratic optimization problem [17,18]. Here, we propose a new method called virtual three-leg supporting algorithm (VTLSA) to obtain the solution of $[\begin{array}{cccc}{F}_{1\mathit{\text{z}}}& F_{2\mathit{\text{z}}}& F_{3\mathit{\text{z}}}& F_{4\mathit{\text{z}}}\end{array}{]}^{\mathrm{T}}$. As shown in Fig.12, the four-legged supporting problem is transformed into four groups of three-legged supporting model: leg 1-2-3 supporting, leg 1-2-4 supporting, leg 1-3-4 supporting, and leg 2-3-4 supporting. Each group is expected to provide a quarter of virtual resultant force in the $z$ direction ($F_{\mathrm{v}\mathrm{r}}$) and virtual resultant torsion in the $x$ and $y$ directions (${\mathit{\tau }}_{\mathrm{v}\mathrm{r}}$). When the virtual $k\text{-}m\text{-}n$ leg supports the lander, we can obtain an analytical solution for the supporting force $[\begin{array}{ccc}{F}_{k\mathit{\text{z}}}^{kmn}& F_{m\mathit{\text{z}}}^{kmn}& F_{n\mathit{\text{z}}}^{kmn}\end{array}{]}^{\mathrm{T}}$ from Eq. (30):

Fig.12 Virtual three-leg supporting algorithm.

Full size|PPT slide

(30)$\left[\begin{array}{ccc}1& 1& 1\\ r_{\mathrm{c}\mathrm{o}\mathrm{m}y}+r_{ky}& r_{\mathrm{c}\mathrm{o}\mathrm{m}y}+r_{my}& r_{\mathrm{c}\mathrm{o}\mathrm{m}y}+r_{ny}\\ r_{\mathrm{c}\mathrm{o}\mathrm{m}x}+r_{kx}& r_{\mathrm{c}\mathrm{o}\mathrm{m}x}+r_{mx}& r_{\mathrm{c}\mathrm{o}\mathrm{m}x}+r_{nx}\end{array}\right]\left[\begin{array}{c}{F}_{k\mathit{\text{z}}}^{kmn}\\ F_{m\mathit{\text{z}}}^{kmn}\\ F_{n\mathit{\text{z}}}^{kmn}\end{array}\right]=\frac{1}{4}\left[\begin{array}{c}{m}_{\mathrm{b}}\left({\ddot{\mathit{\text{z}}}}_{\mathrm{b}}-g\right)\\ I_{\mathrm{b}x}{\dot{\omega }}_{\mathrm{b}x}\\ I_{\mathrm{b}y}{\dot{\omega }}_{\mathrm{b}y}\end{array}\right].$

We can then obtain the expected leg force $F_{i\mathit{\text{z}}}$ by summing the ith vertical force in leg k-m-n supporting ($F_{i\mathit{\text{z}}}^{kmn}$ ($i=k$, m, or n)) for each group, which can be expressed as follows:

(31)$\{\begin{array}{c}{F}_{1\mathit{\text{z}}}=F_{1\mathit{\text{z}}}^{123}+F_{1\mathit{\text{z}}}^{124}+F_{1\mathit{\text{z}}}^{134},\\ F_{2\mathit{\text{z}}}=F_{2\mathit{\text{z}}}^{123}+F_{2\mathit{\text{z}}}^{124}+F_{2\mathit{\text{z}}}^{234},\\ F_{3\mathit{\text{z}}}=F_{3\mathit{\text{z}}}^{123}+F_{3\mathit{\text{z}}}^{134}+F_{3\mathit{\text{z}}}^{234},\\ F_{4\mathit{\text{z}}}=F_{4\mathit{\text{z}}}^{124}+F_{4\mathit{\text{z}}}^{134}+F_{4\mathit{\text{z}}}^{234}.\end{array}$

In accordance with Eq. (14), the real-time adjust tiptoe force ${\mathit{F}}_{\mathrm{t}\mathrm{i}\mathrm{p}}$ can be transformed to the reference joint torque ${\mathit{\tau }}_{\mathrm{r}\mathrm{j}}$.

Currently, the Moore–Penrose inverse matrix can be computed quickly by using matrix decomposition techniques. The method provided by Katsikis and Pappas [33] can be utilized to solve large sparse matrix. The algorithm from Courrieu [34] was based on Cholesky factorization. Toutounian and Ataei [35] presented a CGS-MPi algorithm based on the conjugate Gram–Schmidt process, which is relatively robust for large sparse and rank deficient matrix. Katsikis et al. [36] developed a method based on QR factorization, which functions for singular square matrix or rectangular matrix. An efficient iterative algorithm provided by Esmaeili et al. [37] always converges with fourth-order. Tab.2 lists the computation time of six different algorithms and our method to obtain the tiptoe force with a tolerance of 10⁻¹³. For the Moore–Penrose inverse computation problem of small rectangular matrix that occurs during buffer landings, the methods of generalized inverse (Ginv) [33], tensor product matrix (TPM) [38], and improved Qrginv (IMqrg) [39], have the fastest solution speed. The methods of singular value decomposition (SVD), QR generalized inverse (Qrg) [36], and VTLSA have medium speed, and generalized inverse (Geninv) [34] method has the slowest solution speed. However, the difference among them is small, so they can all be used for stabilizer design of the FLLWR.

Tab.2 Computation time of different algorithms to solve Eq. (29)

Method	Time/s
SVD	3.81 × 10−5
Ginv	1.66 × 10−5
TPM	1.26 × 10−5
Geninv	5.10 × 10−5
Qrg	3.57 × 10−5
IMqrg	2.22 × 10−5
VTLSA	3.34 × 10−5

5 5 Landing numerical simulation

The numerical landing simulations based on the above kinematic and dynamic models are conducted by using the cosimulation of MATLAB and NX to verify the effectiveness of the landing control method and the safety of the driving system. The same parameters, such as the gravitational accelerations on the Moon $g_{\mathrm{m}}=1.633$ m/s², the ground stiffness coefficient K_g = 12000 N/m, and the ground damping coefficient B_g = 2000 N·s/m, are used in different simulations. The buffer landing results on the lunar surface are obviously different under different initial landing conditions: the vertical velocity $v_{\mathit{\text{z}}}$, the horizontal velocity $v_{\mathrm{h}}$, and the load mass. Theoretically, if the lander reaches above the landing site and the landing conditions are satisfied, the robot will finish the landing with a vertical impact. Hence, a vertical landing simulation is designed. When the lander descends from a height of 2.25 m with a total weight of 180 kg, that is, the load weight is 3.5 times of the self-weight, the vertical landing velocity reaches up to $v_{\mathit{\text{z}}}=-2.1\mathrm{m}/\mathrm{s}$ at the touching-ground moment. A parabolic landing simulation ($v_{\mathrm{h}}\ne 0$ m/s) is designed to analyze the landing capacity against a horizontal disturbance velocity. Considering the symmetric distribution of the legs, $v_x=v_{\mathrm{h}}$ is set to reduce the DoF of the lunar gravity testing platform from six to five. The initial conditions in the numerical simulations are listed in Tab.3. These parameters are selected as an example because they are typical conditions. However, other similar parameters can be used.

Tab.3 Initial conditions of numerical simulation

Situation	Weight/kg	vx/(m·s−1)	vz/(m·s−1)
Vertical landing	180	0.0	−2.1
Parabolic landing	180	0.4	−2.1

The results of vertical landing simulation are shown in Fig.13. At $t=1.367\mathrm{s}$, the torques of all legs change intensely, which signifies the occurrence of touching-ground, and the robot velocity in the z direction reaches up to the maximum −2.1 m/s. The legs then begin to compress and absorb the landing impact energy. The torque of the side IDU is about zero because the ground impact energy is mainly distributed in the vertical direction rather than the horizontal direction during the landing. At $t=1.578\mathrm{s}$, the peak torque reaches 155 N·m, and the torque margin is 22.5% relative to the rated peak torque 200 N·m. The thigh IDU torques $T_{\mathrm{t}}$ and the shank IDU torques $T_{\mathrm{s}}$ reach to the stable value after a damping oscillation period. At $t=1.708\mathrm{s}$, $v_{\mathit{\text{z}}}$ is zero, and the robot body reaches the lowest height of 0.3983 m. During the whole landing process, the roll and pitch angles vary in a small range close to zero.

Fig.13 Results of vertical landing simulation: (a) joint torque in each leg, (b) body angle of roll and pitch, and (c) position and velocity of the body.

Full size|PPT slide

Fig.14 illustrates the results of the parabolic landing simulation. At $t=1.367\mathrm{s}$, all legs touch the ground, and the robot states are the same as the ones in the vertical landing simulation, except a horizontal disturbance velocity $v_x=0.4$ m/s. At $t=1.635$ s, the thigh IDU reaches up to peak torque 167.6 $\mathrm{N}\cdot \mathrm{m}$ with a torque margin 16.2%. $T_{\mathrm{t}}$ and $T_{\mathrm{s}}$ reach the stable value after a damping oscillation period. At $t=1.684\mathrm{s}$, the robot body reaches the lowest height of 0.3965 m with $v_z=0$ m/s. Unlike the roll angle that remains in the initial value zero, the pitch angle fluctuates obviously with a maximum value of −7.898°, but it is quickly stabilized to be zero in about 3 s. Compared with the vertical landing simulation, the maximum torque is increased by 12.6 N·m, and the torque margin is decreased by 6.3% correspondingly.

Fig.14 Results of parabolic landing simulation: joint torques in (a) legs 1 and 4 and (b) legs 2 and 3; (c) body angle of roll and pitch; and (d) position and velocity of the body.

Full size|PPT slide

6 6 Prototype and experiments

6.1 6.1 Lunar gravity simulation platform

As shown in Fig.15, a 5-DoF-LGTP, including x, z, $R_x$, $R_y$, and $R_{\mathit{\text{z}}}$ DoFs, is built to simulate the lunar gravity environment. The actual gravity of lander is the same as that on the Moon by using a counterweight system (Fig.16). This platform consists of two pillars, a horizontal guide rail, a cantilever, a wire rope, two fixed pulleys, two counterweights, and sandy soils. A laser tracker is used as the measurement equipment to record the position and velocity variations, and the angle variations are measured by an IMU. The mass of counterweights can be calculated by using the following equations:

Fig.15 Components of the 5-DoF-LGSP.

Full size|PPT slide

Fig.16 Counterweight system design of the 5-DoF-LGSP.

Full size|PPT slide

(32)$\{\begin{array}{c}\left(m_{\mathrm{b}}+m_{\mathrm{c}2}\right)g_{\mathrm{e}}-m_{\mathrm{c}1}{g}_{\mathrm{e}}=\left(m_{\mathrm{b}}+m_{\mathrm{c}1}+m_{\mathrm{c}2}\right)g_{\mathrm{m}},\\ m_{\mathrm{b}}+m_{\mathrm{c}1}+m_{\mathrm{c}2}=m_{\mathrm{t}},\end{array}$

where $m_{\mathrm{b}}=40$ kg is the body mass, $m_{\mathrm{c}1}$ and $m_{\mathrm{c}2}$ are the masses of counterweight 1 and counterweight 2, $m_{\mathrm{t}}$ is the system mass of the lander and loads, and $g_{\mathrm{m}}$ and $g_{\mathrm{e}}$ are the gravitational accelerations on the Moon and Earth, respectively, that is, $g_{\mathrm{m}}=g_{\mathrm{e}}/6$. In this study, we aim to simulate the landing with a 140 kg payload, which indicates $m_{\mathrm{t}}=140$ kg. We obtain $m_{\mathrm{c}1}=75\mathrm{k}\mathrm{g}$ and $m_{\mathrm{c}2}=$ 65 kg.

6.2 6.2 Landing experiments

Similar to numerical simulations, two typical experiments, including the vertical landing and the parabolic landing, are conducted to validate the landing method and analyze the landing performance. The landing initial conditions are the same as those in the simulations.

As a representative example, Fig.17 shows the robot keyframe snapshots in the parabolic landing experiments. (1) is the initial position that is ready to land. (2) is the descending period in the air. (3) is the status at touching ground moment. (4) is compressing to the lowest position (5). (6) is extending to the highest position (7). The robot keeps the stable position after a damping oscillation (8).

Fig.17 Keyframe snapshots in the parabolic landing.

Full size|PPT slide

The vertical landing experiments are conducted, and the results are shown in Fig.18. The lander begins to land from a height of 2.25 m with a weight of 180 kg. At $t=1.313\mathrm{s}$, the torques have a sharp variation, and the lander touches the ground with a vertical velocity $v_{\mathit{\text{z}}}=-2.1\mathrm{m}/\mathrm{s}$. The lander vertical speed is zero, and the body is compressed to the lowest height of 0.3867 m after a deceleration for 0.192 s. The damping oscillation in the z direction lasts about 1 s, and then the body height is stabilized. The peak torque of the thigh IDU reaches up to −166.6 N·m with a margin of 16.7%. The pitch angle has a maximum deviation of −3.135°, but it is quickly adjusted to zero by the attitude stabilizer.

Fig.18 Results of vertical landing experiment: joint torques in (a) leg 1, (b) leg 2, (c) leg 3, and (d) leg 4; (e) body angle of roll and pitch; and (f) position and velocity of the body.

Full size|PPT slide

The parabolic landing experiments are conducted with the same condition as that in vertical landing experiments except for an extra horizontal disturbance velocity. The experiment results are illustrated in Fig.19. At $t=1.46\mathrm{s}$, the torque steps occur, signifying the ground-touching with a horizontal velocity $v_x=0.4\mathrm{m}/\mathrm{s}$ and a vertical velocity $v_{\mathit{\text{z}}}=2.1\mathrm{m}/\mathrm{s}$. The vertical speed becomes zero, and the body is compressed to the lowest height of 0.3512 m after a deceleration for 0.189 s. The damping oscillation in the z direction lasts about 1.142 s, and then the body height is stabilized. The horizontal velocity is stabilized to zero after 1.893 s. The peak torque reaches up to −179.4 N·m with a margin of 10.3%. The roll angle is still about zero, and the pitch angle has a maximum deviation of −8.583°, but it is quickly controlled to be zero by the attitude stabilizer.

Fig.19 Results of parabolic landing experiment: joint torques in (a) leg 1, (b) leg 2, (c) leg 3, and (d) leg 4; (e) body angle of roll and pitch; and (f) position and velocity of the body.

Full size|PPT slide

The above experimental results are consistent with the simulation results. In the vertical landing, the small difference in the maximum torque between the simulation and the experiment is about 11.6 N·m, whereas it is about 11.8 N·m in the parabolic landing. The torque margin is sufficient with minimum value of 10.3% relative to the rated peak torque 200 N·m. In the simulations and experiments, the roll angle changes minimally with only 2° of fluctuation amplitude. The pitch angle changes significantly due to the horizontal impact from the x direction in the parabolic landing, reaching −7.898° and −8.583° in the simulation and experiment, respectively, but it is adjusted to zero in a short time through the height and attitude stabilizer. In the simulation and experiment, the lowest height parameter, which determines whether the engine nozzle at the bottom of the body collides with the ground, is higher than the design requirement (0.2575 m) that is determined by the distance from ground to the bottom of the body (0.2000 m) and the distance from the bottom of the body to the origin of the body (0.0575 m). Specifically, the lowest height in parabolic landing is lower than the one in vertical landing, and the minimum value 0.3512 m in parabolic landing still has a 36.4% safety margin. The proposed landing method is proved to be efficient by the successful landing in the simulations and ground experiments.

7 7 Conclusions and future work

In this study, an FLLWR is designed to execute the landing and roving missions on the Moon. To dissipate impact energy and achieve a safe landing, we propose a landing control method that consists of five parts: landing method framework, state switching strategy, trajectory planning, compliance planning, and height and attitude stabilizer. All sensor data, including encoder, force sensor, and IMU, are fused by the control method framework. The lander can make different responses to events in the landing process via the state switching. The height and attitude stabilizer are designed to realize the expected body height and attitude. A VTLSA is presented to quickly obtain the analytic solution of the stabilizer in the real-time control system. In the dynamic numerical simulation, the lander can achieve stable and safe landing under two typical working conditions: vertical landing and parabolic landing. The proposed control method is proved to be efficient for the landing in the lunar gravity by experiments on the 5-DoF-LGTP with a vertical landing velocity of ‒2.1 m/s and a payload of 140 kg, that is, 3.5 times of the self-weight. With or without an additional horizontal disturbance velocity of 0.4 m/s, the peak torque is 179.4 and 166.6 N·m, respectively, and the height safety margin, which avoids the collision between the engine nozzle installed on the bottom of the robot and the ground, is 36.4% and 50.1%, respectively.

The parameters satisfying the engineering application in compliance control are found by means of numerical dynamics simulation. In the future, we will study the relationship between landing and these parameters, and complete their optimal design systematically. The robot dynamic performance will be further explored to improve the buffer landing.

References

1. Kyung Choi Jin Gil Lee Dai Thermal characteristics ofthe spindle bearing system with a gear located on the bearing spanInternational Journal of Machine Tools & Manufacture 1998 3810171030
2. Bossmanns Bernd Tu Jay F A thermal model for high speedmotorized spindlesInternational Journalof Machine Tools & Manufacture 1999 3913451366
3. Lee Sun-Kyu Yoo Jae-Heung Yang Moon-Su Effect of thermal deformation on machine tool slide guidemotionTribology 2003 364147
4. Yuan J Ni J The real-time error compensationtechnique for CNC machining systemsMechatronics 1998 8359380
5. Yang Hong Ni Jun Adaptive model estimation ofmachine-tool thermal errors based on recursive dynamic modeling strategyInternational Journal of Machine Tools & Manufacture 2005 45111
6. Chen J S Chiou G Quick testing and modelingof thermally induce error of cnc machine toolsInternational Journal of Machine Tools & Manufacture 1995 3510631074
7. Sribastrava A K Veldhuis S C Modelling geometric and thermalerrors in a five-axis CNC machine toolInternationalJournal of Machine Tools & Manufacture 1995 3513211337
8. Chen Zhaonian Chen Zichen The Basic of Machine Tool'sThermal CharacteristicBeijingMechanical Industry Press 1981 (in Chinese)
9. Li Hongqi Yung C Analysis of bearing congurationefects on high speed spindles using an integrated dynamic thermo-mechanicalspindle modelInternational Journal of MachineTools & Manufacture 2004 44347364
10. Fu Jianzhong Chen Zicehn Research on modeling dynamicerrors of precision machine based on fuzzy logic and artifical neuralnetworkJournal of Zhejiang University 2004 63236 (in Chinese)
11. Yang Hong Ni Jun Dynamic neural network modelingfor nonlinear, nonstationary machine tool thermally induced errorInternational Journal of Machine Tools & Manufacture 2005 45455465
12. Yu Zuopin Lu Lu Heat transferingBeijingHighEducation Press 1995 1431 (in Chinese)
13. Zhang Baolin Gu Tongxiang Principles and Methods of NumericalParallel ComputationBeijingNational Defense Industry Press 1999 188189(inChinese)
14. Lapidus Leon Pinder George F Numerical Solution of PartialDifferential Equation in Science and EngineeringSun Nazheng,Translation BeijingCoal Industry Press 1989 3536 (in Chinese)
15. Zhang Chaohui Fan Boqun ANSYS8.0 Thermal AnalysicsBeijingRailroadIndustry Press 2005 (in Chinese)