Landing stability analysis for lunar landers using computer simulation experiments

Abstract

In the overall design and performance prediction of spacecrafts intended to accomplish soft landing on the moon, touchdown dynamics analysis during the landing phase is one of the most important tasks. Due to uncertainties of the landing site and flight states of the vehicle at touchdown, the landing stability is systematically studied in this article by employing the techniques of computer experiments. Firstly, the landing dynamics simulation model of a lunar lander is established and verified by physical tests. Then, design of experiments and response surface methodology have been utilized to analyze the landing stability. Consequently, main effects of initial landing conditions (inputs) on landing dynamics responses (outputs) have been investigated. Response surface models for approximating relationships between inputs and outputs are also acquired. The worst sets of landing conditions and boundaries of initial velocities for stable landings are finally determined with the help of response surface models. The approach for landing stability analysis proposed in this article has been proven effective and efficient and provides valuable guidance for the lander design and analysis.

Keywords

Lunar lander touchdown dynamics landing stability design of experiments response surface methodology

Introduction

Recently, with scientific and technological advancements, many countries or organizations are carrying out extensive research on space and the Moon. The Moon, as the nearest celestial body to the Earth, has become increasingly attractive for space scientists because of its abundant natural resources and high potential as a future territory.¹ Given these benefits of lunar exploration, Americans, Europeans, Japanese, Chinese, and Indians are all planning to go back to the Moon. From the Apollo Program to Lunar CATALYST designed by NASA,² it indicates that scientists and engineers at NASA have never lost interest in lunar exploration.

In lunar exploration, the achievement of soft landing is a challenging task, since the landing operation will occur at an unprepared site, and the vehicle could be oriented randomly with respect to the flight path. Generally, the engine of a lunar lander is shut down at a certain height from the lunar surface to prevent the splash of the lunar regolith, then the lander will freely fall to the lunar surface.³ During the touchdown seconds, great impact forces may be induced to the spacecraft. Therefore, effective measures must be taken to attenuate the shock energy. As reported in the previous lunar exploration missions, such as unmanned Luna, Surveyor, and manned Apollo Program, the legged-type lander is widely used in spacecraft to implement soft landing.⁴ Such landers employ shock absorbers located in the landing gears as the primary shock attenuation component. To assure a successful landing, the landing gear must provide sufficient energy absorption capabilities and also must prevent the landing vehicle from toppling. Thus, in the overall design and performance prediction of the landing gear, touchdown dynamics and landing stability analysis are always of great concern for guaranteeing safe landings of the lunar module and protecting the devices inside it from being damaged.^6,7

Researchers have been dedicated to developing mathematical models to accomplish the soft landing dynamics simulation of legged-type lunar lander in the past decades.^5

–8 The finite element method and the multi-body system dynamics method are two primary modeling approaches. The finite element method is able to consider the nonlinear characteristics of the geometry, materials, and contact so that more accurate results are usually produced. One drawback of this method is that the analysis is too time-consuming. The multi-body dynamics model, despite less precise, can improve the computational efficiency substantially. In this article, a lunar landing dynamics model is built based on the multi-body system dynamics approach and verified by physical tests, which is the foundation for subsequent landing stability analysis. In addition, the analysis of landing stability has received great interest from the literature. Jones and Hinchey⁹ and Black¹⁰ studied the effect of landing condition parameters on stability capability and provided guidelines for the design of landing gears. Alderson and Wells¹¹ and Sperling and Garba¹² created mathematical models and computer programs for investigating the interaction of several parameters on landing stability of Surveyor lunar landers. Hilderman et al.¹³ and Rogers¹⁴ proposed a new analytical method to determine the worst landing conditions and stable landing boundaries for Apollo lunar landers. Drop tests on scaled and full-size landers were also conducted to verify the landing stability capabilities of these landers.^15,16 Landing stability analysis approaches for the Surveyor and the Apollo landers are still in use nowadays.^17

–20 However, the related works mentioned above only considered limited landing conditions and selected values of input factors by arbitrary judgments. No theoretical analysis on the experiment results has been carried out to indicate whether the experiments were sufficient and the consequences were convincing enough. This article proposes a more effective approach by applying design of experiments (DOE) and response surface methodology (RSM) on landing stability analysis. DOE based on the statistical theory has been employed to arrange the landing experiments, which aims to acquire as much information as possible with the least possible effort. Regression analysis is then conducted to infer the influence trend and extent of each input factor on the outputs. The RSM helps provide surrogate models, which have been validated by the analysis of variance (ANOVA), for predicting the extreme value of each response, which contributes to find the worst landing cases more conveniently and convincingly.

The outline of this article is as follows. The first section is the introduction; the landing dynamics model is established and verified in “Landing dynamics model” section; “Computer simulation experiment” section gives the theoretical backgrounds on DOE and RSM; landing stability analysis by using DOE and RSM as well as the results and discussions are finally presented in “Results and discussions” section; the last section is the conclusion.

Landing dynamics model

Configuration of the lunar lander

The lunar lander discussed here is composed of a main body and four landing gear assemblies, which are uniformly distributed around the main body. Each landing gear consists of a primary leg strut, two secondary leg struts, and a footpad, as shown in Figure 1. The leg struts are the primary load-bearing components and serve as supports for the whole lander. Shock absorbers located inside the strut are employed to absorb the landing energy while the lander impacts with the lunar surface. The footpad is beneficial to decreasing the impact penetration on the low bearing strength surface.

Figure 1.

Configuration of the landing gear.

From Figure 1, it can be seen that each leg strut comprises an outer cylinder, an inner cylinder, and a shock absorber. The outer cylinders of all the leg struts are attached to the main body at their upper ends. The inner cylinder of the primary leg strut is connected to the footpad at its lower end, while that of each secondary leg strut is connected to the primary leg strut. Crushable aluminum honeycomb cartridges, which collapse at a constant load, are used as the primary type of shock-attenuation devices.²¹ The primary shock absorber is only capable of absorbing energy during the strut’s compressing, while each secondary shock absorber can absorb energy during the strut’s extending or compressing.

Dynamics model

The landing dynamics simulation model is established by using the multi-body dynamics simulation software ADAMS [version MSC ADAMS 2005]. Each leg strut is attached to the main body through a universal joint with two rotational degrees of freedom. A translational joint is applied to enable the strut’s compressing and extending between the outer and inner cylinders of each leg strut. Each secondary strut is connected to the outer cylinder of the primary strut with a spherical joint permitting all three rotations, and the joint connecting the footpad and the inner cylinder of the primary strut is also a spherical type.

The significant forces acting on the landing vehicle include the gravity, the shock absorber forces, and the contact forces between the footpads and the lunar surface. The honeycomb shock absorber in each strut is assumed to produce acting and reacting forces directing along the axis of the strut at the outer and inner cylinders of the strut. Figure 2 illustrates typical honeycomb configurations for primary and secondary strut assemblies. Note that two stages of load level are designed for the primary strut. The initial low force level is intended to enhance the toppling stability, while the higher force level limits loads on the vehicle body without bottoming out the strut.²²

Figure 2.

Honeycomb configurations for primary and secondary struts. Load-deflection curves of honeycombs in the (a) primary strut and (b) secondary strut.

The shock absorber force of the primary strut F _P is mathematically expressed as

F_{P} = {\begin{matrix} F_{HP} (S_{P}) - c_{v P} \cdot {\dot{S}}_{P}, {\dot{S}}_{P} > 0 \\ 0, {\dot{S}}_{P} \leq 0 \end{matrix}

where F _HP(S _P) is the honeycomb force as the primary stroke varies, whose value equals F _HP1 when 0 < S _P ≤ S _HP1 and takes F _HP2 when S _HP1 < S _P ≤ S _HP2, as illustrated in Figure 2(a). c_v _P represents the vicious damping coefficient of the primary strut, and ${\dot{S}}_{P}$ is the velocity of the inner cylinder relative to the outer cylinder of the primary strut.

The shock absorber force of each secondary strut F _S is computed by

F_{S} = F_{HS} (S_{S}) - c_{v S} \cdot {\dot{S}}_{S}

where F _HS(S _S) represents the honeycomb force as the secondary stroke varies, whose value equals F _HSC when 0 < S _S ≤ S _HSC and takes F _HSE when S _HSE ≤ S _S < 0, as illustrated in Figure 2(b). c_v _S represents the vicious damping coefficient of the secondary strut and ${\dot{S}}_{S}$ is the relative velocity between the inner and outer cylinders of the secondary strut.

The contact force between the footpad and the lunar surface is divided into a normal force perpendicular to the lunar surface and a tangential force parallel to the surface. The normal force F_n is modeled as a nonlinear spring damper²³

F_{n} = K \cdot δ^{e} + D \cdot \dot{δ}

in which δ is the penetration of the footpad into the lunar surface and $\dot{δ}$ is the derivative of the penetration. K represents the contact stiffness coefficient, e represents the nonlinear exponent, and D represents the damping coefficient. The stiffness component models the elasticity of the contact surface, and the damping component reflects the dissipation of energy. The tangential force F_t is mathematically modeled as follows

F_{t} = μ \cdot F_{n} + f_{v} \cdot v_{t}

where μ is the Coulomb friction coefficient, v_t is the tangential velocity of the footpad, and f_v ⋅ v_t represents the friction force due to velocity.

Physical landing tests have been conducted on a 1/4 lander with only a single leg. The corresponding dynamics model is also built in ADAMS (Figure 3). Results of tests are compared with those of simulations to verify the effectiveness of the modeling approach mentioned above. Physical tests are implemented under two cases with different landing weights and the same release height as shown in Table 1. Figures 4 and 5 illustrate the comparisons of the time history curves of the primary and secondary strokes between tests and simulations under the two landing cases. It can be seen that simulation results are in good agreement with test results, which validates the correctness of our dynamics model. The single-leg model is then expanded to the complete four-leg model (Figure 6), ready for use in subsequent landing stability analysis.

Figure 3.

Dynamics simulation model for the single-leg lander.

Table 1.

Cases for landing tests of the single-leg lander.

Case no.	Release altitude (m)	Landing weight (kg)
S1	1.0	330
S2	1.0	490

Figure 4.

Comparison of test and simulation results under case S1. Time history curves of the (a) primary stroke and (b) secondary stroke.

Figure 5.

Comparison of test and simulation results under case S2. Time history curves of the (a) primary stroke and (b) secondary stroke.

Figure 6.

Dynamics simulation model for the four-leg lander.

Physical landing tests on the complete four-leg model are also carried out under eight different cases as listed in Table 2. Figure 7 illustrates a schematic diagram for physical landing tests of the four-leg lander.²⁴ Landing case parameters include the release height, the initial horizontal velocity, the slope of the landing surface, and the landing manner. The release height is the vertical distance from the landing surface to the bottom of the footpad that lands first. The initial horizontal velocity is directed toward the slope. A positive value denotes a downhill landing, while a negative value denotes an uphill landing. Two symmetrical landing manners are considered in the tests, known as the 1-2-1 and the 2-2 landing manner. The former represents initial touchdown of one leg, followed by simultaneous touchdown of the middle two legs, and then subsequent touchdown of the remaining legs, as shown in Figure 8(a). The latter represents initial simultaneous touchdown of two legs, followed by simultaneous touchdown of the other two legs, as shown in Figure 8(b).

Table 2.

Cases for landing tests of the four-leg lander.

Case No.	Release height (m)	Horizontal velocity (m s⁻¹)	Slope (°)	Landing mode
C1	0.8	0.0	0	1-2-1
C2	1.0	1.0	0	1-2-1
C3	1.0	0.0	8	1-2-1
C4	1.0	0.0	8	2-2
C5	1.0	1.0	8	1-2-1
C6	1.0	1.0	8	2-2
C7	1.0	−1.0	8	1-2-1
C8	1.0	−1.0	8	2-2

Figure 7.

Schematic diagram for physical landing tests of the four-leg lander.

Figure 8.

Schematic diagrams for symmetrical landing manners. (a) 1-2-1 landing manner and (b) 2-2 landing manner.

According to the landing cases listed in Table 2, dynamics simulations are also executed with the model shown in Figure 6. Results of the maximum primary strokes for all legs under each case from simulations are listed in Table 3 and compared with those from tests. It can be seen that most of the simulation results are close to the corresponding test results. Significant differences in several cases are mainly caused by three major factors. Firstly, the dynamics simulation model is incapable of being completely consistent with the actual scenario. In particular, it is very difficult for accurately simulating the nonlinear forces between footpads and the landing surface. Secondly, fabrication deviations may cause test results different from theoretical ones. Especially, since honeycomb forces play great roles in the landing process, manufacturing precisions of the honeycombs (which are nonreusable and diverse honeycombs have been used in different test cases) may have a big impact on test results. Thirdly, the inaccuracy of observation or measurement is another source of error for physical tests. In conclusion, we consider that the simulation model is accurate enough for prediction of landing dynamics responses.

Table 3.

Comparisons of maximum primary stroke results between simulations and tests.

Case no.		Leg 1 SPC₁ (mm)	Leg 2 SPC₂ (mm)	Leg 3 SPC₃ (mm)	Leg 4 SPC₄ (mm)
C1	Simulation	121	121	121	121
C1	Test	119	114	104	128
C2	Simulation	105	151	210	151
C2	Test	113	154	207	161
C3	Simulation	183	123	116	123
C3	Test	171	138	111	144
C4	Simulation	164	90	90	164
C4	Test	160	118	119	150
C5	Simulation	240	133	77	133
C5	Test	240	136	83	123
C6	Simulation	192	77	77	192
C6	Test	201	69	98	207
C7	Simulation	142	94	146	94
C7	Test	144	83	155	107
C8	Simulation	134	110	110	134
C8	Test	153	111	126	139

Computer simulation experiment

To analyze the landing stability, the techniques of computer simulation experiment are introduced in this section, which are composed of the DOE and the RSM.

Design of experiments

In general, DOE is a statistical method for arranging experiments and analyzing experimental data. It is an effective strategy for examining the behavior of a physical phenomenon when changing its input values in order to detect their effects on the output.²⁵ DOE helps in understanding the input–output relationship and also in identifying the relative importance of the inputs.

To start with, some fundamental terms involved in the DOE technique are introduced. The factors are the inputs of an experiment, which are varied during the DOE. A factor may be assigned at least two distinct values in the DOE, and each value is called a level. In other words, levels are the different possible values taken from each factor. The response is the output of an experiment, which reflects how the system responds to a certain combination with different values of input factors. The relationship between the response y and the factor vector x = (x ₁,…, x_r ) ^T can be written by

y = y (x)

Practically, the mean of the response is calculated, after the implementations of n experiments, as follows

E (y) \approx \bar{y} (D_{n}) = \frac{1}{n} \sum_{i = 1}^{n} y (x_{i})

where D _n is the design matrix, and each row represents a sample point in an r-dimensional space as shown below

D_{n} = [\begin{matrix} x_{11} & x_{12} & \dots & x_{1 r} \\ x_{21} & x_{22} & \dots & x_{2 r} \\ \cdot & \cdot & \dots & \cdot \\ \cdot & \cdot & \dots & \cdot \\ \cdot & \cdot & \dots & \cdot \\ x_{n 1} & x_{n 2} & \dots & x_{n r} \end{matrix}]

In DOE, the most important task is how to construct the design matrix D _n . The optimal Latin hypercube sampling (LHS) is one of the most widely used sampling approach in the design of computer experiments, which enhances the space-filling property of the sample points.²⁶ The optimal LHS employs the max–min distance for maximizing the minimum distance among all distances between sample points. The optimal criterion can be represented as

min_{1 \leq i, j \leq n, i \neq j} d (x_{i}, x_{j}) \to max

where d( x _i , x _j ) represents the Euclidean distance between two sample points x _i and x _j , defined as

d (x_{i}, x_{j}) = d_{i j} = {[\sum_{m = 1}^{k} {| x_{i m} - x_{j m} |}^{t}]}^{1 / t}, t = 1 or 2

Response surface methodology

The RSM provides an approximate mathematical expression for the relationship between a response y and its factors x ₁, x ₂,…, x_r through regression analysis.²⁷ Consider the following second-degree polynomial expression as the approximation model

y = β_{0} + \sum_{i = 1}^{r} β_{i} x_{i} + \sum_{i = 1}^{r} β_{i i} x_{i}^{2} + \sum_{i j} \sum_{(i < j)} β_{i j} x_{i} x_{j} + ε

The model contains (r + 1)(r + 2)/2 regression parameters, including coefficients for main effects β_i , coefficients for quadratic main effects β_ii , and coefficients for interaction effects β_ij . ε is the fitting error and assumed to follow a normal distribution. In order to estimate the unknown coefficients of equation (10), a series of experiments should be carried out, and the response is observed after each experiment. This creates a set of equations in the following form

y = U β + ε

where y = [y ₁, y ₂,…, y_n ] ^T is the response vector in n experiments and ε = [ε ₁, ε ₂,…, ε_n ] ^T is the error vector. β is the coefficient vector and given by

β = {[\begin{matrix} β_{0} & β_{1} & \dots & β_{r} & β_{11} & \dots & β_{r r} & β_{12} & \dots & β_{r - 1 r} \end{matrix}]}^{T}

X represents the input vector as follows

X = {[\begin{matrix} 1 & x_{1} & \dots & x_{r} & x_{1}^{2} & \dots & x_{r}^{2} & x_{1} x_{2} & \dots & x_{r - 1} x_{r} \end{matrix}]}^{T}

U is the design matrix having the following form

U = [\begin{matrix} 1 & x_{11} & \dots & x_{1}_{r} & x_{11}^{2} & \dots & x_{1 r}^{2} & x_{11} x_{12} & \dots & x_{1 r - 1} x_{1 r} \\ 1 & x_{21} & \dots & x_{2}_{r} & x_{21}^{2} & \dots & x_{2 r}^{2} & x_{21} x_{22} & \dots & x_{2 r - 1} x_{2 r} \\ . & \dots & \dots & \dots \\ . & \dots & \dots & \dots \\ . & \dots & \dots & \dots \\ . & \dots & \dots & \dots \\ 1 & x_{n 1} & \dots & x_{n}_{r} & x_{n 1}^{2} & \dots & x_{n r}^{2} & x_{n 1} x_{n 2} & \dots & x_{n r - 1} x_{n r} \end{matrix}]

Using the least squares regression, the estimate of β can be obtained as follows

\hat{β} = {(U^{T} U)}^{- 1} U^{T} y

Furthermore, each component of β is estimated as follows

{\begin{matrix} {\hat{β}}_{0} = \frac{1}{n} \sum_{u = 1}^{n} y_{u} \\ {\hat{β}}_{i} = \sum_{u = 1}^{n} x_{u i} y_{u} / \sum_{u = 1}^{n} x_{u i}^{2} \\ {\hat{β}}_{i i} = \sum_{u = 1}^{n} x_{u i}^{2} y_{u} / \sum_{u = 1}^{n} {(x_{u i}^{2})}^{2} \\ {\hat{β}}_{i j} = \sum_{u = 1}^{n} x_{u i} x_{u j} y_{u} / \sum_{u = 1}^{n} {(x_{u i} x_{u j})}^{2} \end{matrix}

By differentiating equation (10), the first-, second-order main effects, and interaction effects can be calculated as follows

M_{i} = β_{i} d x_{i}

M_{i i} = 2 β_{i i} {\bar{x}}_{i} d x_{i}

M_{i j} = β_{i j} ({\bar{x}}_{j} d x_{i} + {\bar{x}}_{i} d x_{j})

After obtaining the response surface model, ANOVA is performed to test the significance and the goodness of fit of the model.²⁸ In ANOVA, F-test is usually employed to test the significance of regression. For significance test, we suppose the hypothesis

H_{0} : β_{1} = \dots = β_{r} = β_{11} = \dots = β_{r r} = β_{12} = \dots = β_{r - 1 r} = 0

If H ₀ is accepted, it indicates that the regression is not significant, which means that it is not appropriate to use model (10) to formulate the relationship between the input vector X and the response y. To define the statistic for H ₀ test, we first introduce the total sum of squares of real values from the overall mean

SST = \sum_{i = 1}^{n} {(y_{i} - \bar{y})}^{2}

where y_i (i = 1, 2,…, n) are real observed values of the response in the n experiments, and $\bar{y} = \frac{1}{n} \sum_{i = 1}^{n} y_{i}$ is the overall mean of the n observed data. SST can be further expressed as the sum of regression sum of squares SSR and residual sum of squares SSE

SST = SSR + SSE = \sum_{i = 1}^{n} {({\hat{y}}_{i} - \bar{y})}^{2} + \sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}

in which ${\hat{y}}_{i}$ (i = 1, 2,…, n) are the n predicted values computed by the regression model. SSR represents the variation attributed to the regression, that is, the relationship between the input vector X and the response y. SSE represents the variation due to the error. Now, the statistic for testing H ₀ can be defined as

F = \frac{SSR / m}{SSE / (n - m - 1)}

where m is the number of regression coefficients excluding the constant term and called SSR’s degrees of freedom. (n − m − 1) is called SSE’s degrees of freedom. It can be proved that F follows an F-distribution with parameters m and (n − m − 1), that is, F ∼ F(m, n − m − 1). We then compare the ratio F with the 100(1 − α)% point of the tabulated F(m, n − m − 1) distribution. If F > F_α (m, n − m − 1) (α is termed as the significance level), we reject the hypothesis H ₀ and infer that the regression model is significant running a risk of less than 100α% of being wrong.

A useful statistic for checking how well the regression model approximates the real observed data is the coefficient of determination R ², which is defined as

R^{2} = \frac{SSR}{SST}

R ² measures the proportion of total variation about the mean explained by the regression. The closer R ² value is to 1, the better fit of the model will be. An R ² of 1 indicates that the regression model perfectly fits the data.

Results and discussions

In “Landing dynamics model” section and “Computer simulation experiment” section, we have established the dynamics model of the lunar lander and introduced the theoretical background on computer simulation experiments. In this section, the techniques of DOE and RSM are applied for analyzing the landing stability of the lunar lander. Results and discussions are then presented to demonstrate the effectiveness of the proposed method. Firstly, criteria for judging a safe landing of the lunar lander are specified. After that the main effect analysis is conducted to investigate the trend and the extent of each input factor on the landing stability responses. Then, the worst landing condition for each response is acquired by using the RSM. The envelope of initial velocities for stable landing is given at last.

Landing stability criteria

According to the requirements of lander toppling prevention and impact energy absorption, criteria for stable landing are given as follows:

no vehicle toppling happens;

no excessive stroke of any primary shock absorber occurs; and

no excessive compression or extension stroke of any secondary shock absorber takes place.

A parameter representation of lander toppling is defined as the minimum stability distance D _min. Stability distances are measured from the mass center of the lander to overturn planes, which are vertical planes that pass through every two adjacent landing gear footpads as shown in Figure 9. The stroke of each primary shock absorber is denoted as SPC _i (i = 1, 2, 3, 4), and SPC_max denotes the maximum value taken from results of the four primary strokes. The compression stroke and the extension stroke of the secondary shock absorber are denoted as SSC _i and SSE _i (i = 1, 2, 3, 4), respectively. SSC_max and SSE_max denote the maximum absolute values among results of the four secondary compression and the four extension strokes, respectively

Figure 9.

Parameter for measuring lander toppling.

Hence, the landing stability criteria can be also expressed as follows:

D _min > 0, where D _min = min(D ₁, D ₂, D ₃, D ₄);

SPC_max ≤ SPC_th, where SPC_max = max(SPC₁, SPC₂, SPC₃, SPC₄);

SSC_max ≤ SSC_th, where SSC_max = max(SSC₁, SSC₂, SSC₃, SSC₄) and SSE_max ≤ SSE_th, where SSE_max = max(SSE₁, SSE₂, SSE₃, SSE₄).

SPC_th, SSC_th, and SSE_th denote threshold values for the stroke of the primary shock absorber, the compression stroke of the secondary shock absorber, and the extension stroke of the secondary shock absorber, respectively.

Main effect analysis

Concerning the landing stability analysis, input factors are selected as follows: the initial vertical velocity v_v , the initial horizontal velocity v_h , the initial yaw angle θ_x , the initial pitch angle θ_y , the lunar slope angle α_s , and the landing weight m_L , as illustrated in Figure 10. Output responses include the minimum toppling distance D _min, the maximum primary compressing stroke SPC_max, and the maximum secondary compression and extension stokes, SSC_max and SSE_max.

Figure 10.

Input factors for landing stability analysis. (a) Front view and (b) top view.

Firstly, the main effect analysis is conducted to investigate the trend and the extent of each input factor on the landing stability responses. Table 4 shows the lower and upper limits of all input factors, which are determined by the guidance and control systems. A series of landing dynamics simulation experiments have been designed by employing the optimal LHS technique introduced in “Design of experiments” section. With the results of simulation experiments and by utilizing the regression analysis of equations (15) to (18), main effects of input factors on each landing stability response are finally obtained and demonstrated in Figures 11 to 14.

Table 4.

Ranges of input factors for main effects analysis.

v_v (m·s⁻¹)	v_h (m·s⁻¹)	θ_x (°)	θ_y (°)	α_S (°)	m_L (kg)
[1, 4]	[−0.8, 0.8]	[0, 45]	[−4, 4]	[0, 10]	[1150, 1300]

Figure 11.

Main effect of each factor on D _min. (a) v_v on D _min, (b) v_h on D _min, (c) θ_x on D _min, (d) θ_v on D _min, (e) α_s on D _min, and (f) m_L on D _min.

Figure 12.

Main effect of each parameter on SPC_max. (a) v_v on SPC_max, (b) v_h on SPC_max, (c) θ_x on SPC_max, (d) θ_v on SPC_max, (e) α_s on SPC_max, and (f) m_L on SPC_max.

Figure 13.

Main effect of each parameter on SSC_max. (a) v_v on SSC_max, (b) v_h on SSC_max, (c) θ_x on SSC_max, (d) θ_v on SSC_max, (e) α_s on SSC_max, and (f) m_L on SSC_max.

Figure 14.

Main effect of each parameter on SSE_max. (a) v_v on SSE_max, (b) v_h on SSE_max, (c) θ_x on SSE_max, (d) θ_v on SSE_max, (e) α_s on SSE_max, and (f) m_L on SSE_max.

Figure 11 shows the main effect of each factor on the response D _min. It appears that in an unsymmetrical downhill landing manner (0° < θ_x < 45° and v_h > 0) and with the minimum landing weight, the lander is more likely to overturn as v_v , v_h , θ_v , and α_s keep increasing.

Figure 12 shows the main effect of each factor on the response SPC_max. It can be seen that in a 1-2-1 uphill landing manner (θ_x = 0° and v_h < 0) and with the maximum landing weight, the demand for energy absorption of the primary strut will become more critical as the absolute values of v_v , v_h , θ_v , and α_s keep increasing.

Figures 13 and 14 show the main effects of each factor on the responses SSC_max and SSE_max, respectively. It seems that in an unsymmetrical uphill landing manner and with the maximum landing weight, situations for the secondary struts absorbing energy will become worse as the absolute values of v_v , v_h , θ_v , and α_s keep increasing.

Worst landing conditions

Landing dynamics simulations are so complicated and time-consuming that it is extremely difficult to find the worst landing conditions from the infinite combinations of different values of input factors. Therefore, we employ the RSM to approximate the relationships between landing dynamics responses and input variables. Since the approximation model is always more accurate if ranges of the input variables are set small enough, we further narrow the range of each input factor based on the results of the main effect analysis, as listed in Table 5.

Table 5.

Ranges of input factors for response surface models.

Factors Responses	v_v (m·s⁻¹)	v_h (m·s⁻¹)	θ_x (°)	θ_y (°)	α_S (°)	m_L (kg)	Footpad state
D _min	[2, 4]	[0.4, 0.8]	[15, 30]	[2, 4]	[6, 10]	[1150, 1250]	no slip
SPC_max	[2, 4]	[−0.8, −0.4]	[0, 12]	[2, 4]	[7, 10]	[1200, 1300]	no slip
SSC_max	[2, 4]	[−0.8, −0.4]	[9, 24]	[2, 4]	[7, 10]	[1200, 1300]	no slip
SSE_max	[2, 4]	[−0.8, −0.4]	[0, 16]	[2, 4]	[4, 8]	[1200, 1300]	slip

For each response, 128 simulation experiments are arranged by using the optimal LHS approach. Levels of each input factor are all located in the ranges listed in Table 5. Under every single sample point with a combination of specified values of input factors, the landing dynamics simulation is performed. Four sets of simulation results are eventually generated, which are then used to determine the unknown regression coefficients of the four response surface models by applying the least square method, respectively. Consequently, the second-order polynomial response surface model for each landing dynamics response is obtained, as shown in the following equations

\begin{array}{l} D_{min} = 218.408 - 104 .102 v_{v} + 171.573 v_{h} - 14.795 θ_{x} \\ - 179.803 θ_{y} + 1.099 α_{S} + 2.731 m_{L} + 5.976 v_{v}^{2} \\ - 95.249 v_{h}^{2} + 0.553 θ_{x}^{2} + 2.346 θ_{y}^{2} - 1.109 α_{S}^{2} \\ - 0.0015 m_{L}^{2} - 126.149 v_{v} v_{h} - 0.216 v_{v} θ_{x} \\ - 18.487 v_{v} θ_{y} - 31.341 v_{v} α_{S} + 0.258 v_{v} m_{L} \\ - 5.732 v_{h} θ_{x} - 3.183 v_{h} θ_{y} - 25.42 v_{h} α_{S} \\ + 0.397 v_{h} m_{L} + 0.853 θ_{x} θ_{y} + 0.185 θ_{x} α_{S} \\ - 0.0085 θ_{x} m_{L} + 4.833 θ_{y} α_{S} + 0.114 θ_{y} m_{L} \\ + 0.024 α_{S} m_{L} \end{array}

{SPC}_{max} = 765.272 + 58.674 v_{v} + 153.503 v_{h} - 2.071 θ_{x} - 3.716 θ_{y} - 5.98 α_{S} - 1.322 m_{L} - 5.14 v_{v}^{2} + 14.63 v_{h}^{2} - 0.0094 θ_{x}^{2} + 0.156 θ_{y}^{2} - 0.129 α_{S}^{2} + 0.000464 m_{L}^{2} - 6.261 v_{v} v_{h} + 0.065 v_{v} θ_{x} + 1.786 v_{v} θ_{y} + 0.482 v_{v} α_{S} + 0.023 v_{v} m_{L} + 0.056 v_{h} θ_{x} + 4.097 v_{h} θ_{y} - 1.262 v_{h} α_{S} - 0.148 v_{h} m_{L} - 0.008 θ_{x} θ_{y} + 0.0052 θ_{x} α_{S} + 0.0015 θ_{x} m_{L} - 0.227 θ_{y} α_{S} + 0.0038 θ_{y} m_{L} + 0.0053 α_{S} m_{L}

{SSC}_{max} = 467.328 - 77.746 v_{v} + 165.34 v_{h} + 0.363 θ_{x} - 4.447 θ_{y} - 4.189 α_{S} - 0.505 m_{L} + 4.612 v_{v}^{2} + 40.044 v_{h}^{2} - 0.0093 θ_{x}^{2} + 0.145 θ_{y}^{2} - 0.011 α_{S}^{2} + 0.000144 m_{L}^{2} - 9.128 v_{v} v_{h} + 0.09 v_{v} θ_{x} + 1.619 v_{v} θ_{y} + 0.806 v_{v} α_{S} + 0.034 v_{v} m_{L} + 0.065 v_{h} θ_{x} - 3.797 v_{h} θ_{y} + 1.484 v_{h} α_{S} - 0.09 v_{h} m_{L} - 0.0525 θ_{x} θ_{y} - 0.094 θ_{x} α_{S} + 0.000291 θ_{x} m_{L} - 0.145 θ_{y} α_{S} + 0.00076 θ_{y} m_{L} + 0.0048 α_{S} m_{L}

{SSE}_{max} = 131.518 - 25.553 v_{v} + 21.766 v_{h} + 0.563 θ_{x} - 0.538 θ_{y} - 0.28 α_{S} - 0.172 m_{L} + 3 v_{v}^{2} + 13.65 v_{h}^{2} - 0.0028 θ_{x}^{2} + 0.069 θ_{y}^{2} - 0.0213 α_{S}^{2} + 0.0000558 m_{L}^{2} - 7.702 v_{v} v_{h} - 0.196 v_{v} θ_{x} + 0.929 v_{v} θ_{y} + 0.1 v_{v} α_{S} + 0.02 v_{v} m_{L} + 0.8 v_{h} θ_{x} + 0.806 v_{h} θ_{y} + 1.707 v_{h} α - 0.02 v_{h} m_{L}_{S} - 0.018 θ_{x} θ_{y} + 0.03 θ_{x} α_{S} + 0.000042 θ_{x} m_{L} - 0.2215 θ_{y} α_{S} - 0.00096 θ_{y} m_{L} + 0.000025 α_{S} m_{L}

Subsequently, ANOVA is implemented for each regression model, and the results are shown in Tables 6 to 9 (SS represents sum of squares and df represents degrees of freedom). Given the significance level α = 0.01, from the F-distribution table, we have the critical value F _0.01(27, 100) = 1.93, which is much smaller than the F-value of each regression model. This means we have 99% confidence that all regression models are significant. It can also be seen from the tables that all R ² values are highly close to 1, which indicates each response surface model fits the relationship between its output and inputs well.

Table 6.

ANOVA table for RSM of D _min.

Variation source	SS	df	F-value	R ²
Regression	3.186E6	27	1043.3	0.996
Residual	1.131E4	100	1043.3	0.996
Total	3.187E6	127

RSM: response surface methodology; ANOVA: analysis of variance; SS: sum of squares; df: degrees of freedom.

Table 7.

ANOVA table for RSM of SPC_max.

Variation source	SS	df	F-value	R ²
Regression	2.198E5	27	721.762	0.994
Residual	1127.9	100	721.762	0.994
Total	2.2093E5	127

RSM: response surface methodology; ANOVA: analysis of variance; SS: sum of squares; df: degrees of freedom.

Table 8.

ANOVA table for RSM of SSC_max.

Variation source	SS	df	F-value	R ²
Regression	7776	27	82.797	0.946
Residual	347.84	100	82.797	0.946
Total	8123.84	127

RSM: response surface methodology; ANOVA: analysis of variance; SS: sum of squares; df: degrees of freedom.

Table 9.

ANOVA table for RSM of SSE_max.

Variation source	SS	df	F-value	R ²
Regression	27,404	27	2785.5	0.998
Residual	36.437	100	2785.5	0.998
Total	27,440.437	127

RSM: response surface methodology; ANOVA: analysis of variance; SS: sum of squares; df: degrees of freedom.

Taking advantage of the response surface models, the worst landing condition for each response can then be obtained by finding extremal solutions of equations (25) to (28) via a certain optimization algorithm, such as the sequential quadratic programming algorithm. Figure 15 demonstrates the results of the worst landing conditions for all landing dynamics responses, and Table 10 lists both predicted and real values of each response under the corresponding worst landing case. It seems that the deviation of each predicted value from its real value is acceptable, which proves the feasibility and effectiveness of the RSM.

Figure 15.

Worst sets of landing conditions. Worst landing condition for (a) vehicle toppling, (b) primary stroke, (c) secondary compression stroke, and (d) secondary extension stroke.

Table 10.

Landing dynamic response under the worst condition.

	D _min (mm)	SPC_max (mm)	SSC_max (mm)	SSE_max (mm)
Predicted value	424	203	71	86
Actual value	343	195	83	82
Requirements	>0	<350	<100	<120

Landing stability boundaries

Initial velocities have great influences on the landing stability and are also the essential factors for determining when to shut down the engine and to fall onto the lunar ground. Hence, we attempt to find the envelope of initial velocities for assuring stable landings. Under each worst landing condition, we gradually increase the initial vertical velocity (v_v ) and horizontal velocity (v_h ) simultaneously while keeping other parameters constant. We keep performing the landing dynamics simulation until the lander overturns or any of the strokes exceeds the threshold value. Finally, the boundary curve for each landing response is acquired and plotted with a dotted line as shown in Figure 16. Subject to the constraints of the four boundary curves, the resulting overall boundary for stable landings is then constructed and plotted with solid lines as shown in Figure 16.

Figure 16.

Landing stability boundaries.

The results show that, to ensure a stable landing in the premise that α_S ≤ 10° and θ_y ≤ 4°, values of the initial vertical and horizontal velocities should be within the pentagon envelope in Figure 16. The maximum absolute values of the initial vertical and horizontal velocities should never exceed 5.6 and 1.1 m s⁻¹, respectively.

Conclusion

In this article, an approach for studying the landing stability of lunar landers by employing the techniques of DOE and RSM has been proposed. Some useful conclusions are drawn as follows:

A landing dynamics simulation model for the lunar lander has been built and verified by physical tests. The model enables us to evaluate the landing performances of the lander conveniently and is capable of performing plenty of landing experiments that are difficult or expensive to implement in reality.

Based on the dynamics model, DOE is employed to arrange landing simulation experiments rationally. Trends and extents of input factors, such as v_v , v_h , θ_x , θ_y , α_s , and m_L , on output responses D _min, SPC_max, SSC_max, and SSE_max have been investigated with the least possible effort. The RSM has been proven feasible and effective to approximate the relationships between input factors and output landing dynamics responses. The worst landing condition for each landing stability response has been found with the help of response surface models.

Under the worst landing conditions, landing stability boundaries measured by D _min, SPC_max, SSC_max, and SSE_max have been obtained, subject to which the envelope of initial vertical and horizontal velocities for stable landings is finally constructed. As for the lunar lander we study in this article, to assure stable landings, the initial vertical and horizontal velocities must not exceed 5.6 and 1.1 m s⁻¹, respectively, in the premise that the slope α_S and the pitch angle θ_y are less than 10° and 4°, respectively.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflict of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (11202010 and 61340044).

References

Pham

Zhao

Goo

. Landing stability simulation of a 1/6 lunar module with aluminum honeycomb dampers. Int J Aeronaut Space 2013; 14(4): 356–368.

Shao

Chen

. Trajectory optimization for lunar soft landing with a Hamiltonian-based adaptive mesh refinement strategy. Adv Eng Soft 2016; 100: 266–276.

Chu

Zhang

. Lunar probe techniques. Beijing: China Science and Technology Press [in Chinese], 2007.

Liu

Huang

. Review of the development of soft-landing mechanisms for lunar explorations. Chine Space Sci Technol [in Chinese] 2006; 26(1): 33–39.

Hilchenbach

Kuchemann

Rosenbauer

. Impact on a comet: Rosetta lander simulations. Planet Space Sci 2000; 48(5): 361–369.

Nohmi

Miyahara

. Modeling for lunar lander by mechanical dynamics software. In: AIAA modeling and simulation technologies conference and exhibit, San Francisco, CA, United States, 15–18 August 2005. AIAA. DOI: 10.2514/6.2005-6416.

Wan

Nie

Chen

. Modeling and simulation of lunar lander soft-landing using transient dynamics approach. In: Proceedings of 2010 international conference on computational and information sciences, Chengdu, China, 17–19 December 2010, pp. 741–744. IEEE.

Liang

Chai

HY,

Chen

. Review of soft landing dynamics modeling and analysis for lunar landers. Space Eng [in Chinese] 2011; 20(6): 104–112.

Jones

Hinchey

. Some basic guidelines for establishing structural design parameters for the landing gear of stable, soft landing spacecraft. In: 9th structural dynamics and materials conference, Palm Springs, CA, United States, 1–3 April 1968. AIAA. DOI: 10.2514/6.1968-345.

10.

Black

. Quadrupedal landing gear systems for spacecraft. J Space Rocket 1964; 1(2): 196–203.

11.

Alderson

Wells

. Final report on surveyor lunar touchdown stability study. NASA Technical Report NASA-CR-80514, Pasadena, CA, United States: Jet Propulsion Laboratory, California Institute of Technology, 8 July 1966.

12.

Sperling

Garba

. A treatise on the Surveyor lunar landing dynamics and an evaluation of pertinent telemetry data returned by Surveyor I. NASA Technical Report NASA-CR-87421, Pasadena, CA, United States: Jet Propulsion Laboratory, California Institute of Technology, 15 August 1967.

13.

Hilderman

Mueller

Mantus

. Landing dynamics of the lunar excursion module. J Space Rocket 1966; 3(10): 1484–1489.

14.

Rogers

. Apollo experience report: lunar module landing gear subsystem. NASA Technical Report NASA-TN-D-6850, Houston, TX, United States: NASA Johnson Space Center, 1 June 1972.

15.

Herr

Leonard

. Dynamic model investigation of touchdown stability of lunar-landing vehicles. NASA Technical Report NASA-TN-D-4215, Hampton, VA, United States: NASA Langley Research Center, 1 October 1967.

16.

Blanchard

. Full-scale dynamic landing-impact investigation of a prototype lunar module landing gear. NASA Technical Report NASA-TN-D-5029, Hampton, VA, United States: NASA Langley Research Center, 1 March 1969.

17.

Zhu

Yang

. Touchdown stability simulation of landing gear system for lunar lander. J Astro [in Chinese] 2009; 30(5): 1792–1796.

18.

Briere

Weber

. Robotic lunar lander development project: three dimensional dynamic stability testing and analysis. In: Proceedings of 26th space simulation conference, Annapolis, MD, United States, 18–21 October 2010. NASA/CP-2011-215873.

19.

Ponnusamy

Maahs

. Development and testing of leg assemblies for robotic lunar lander. In: Proceedings of 14th European space mechanisms and tribology Symposium, Constance, Germany, 28–30 September 2011, pp. 485–492. ESA communications.

20.

Cole

Bassler

Cooper

. The challenges of designing a lightweight spacecraft structure for landing on the lunar surface. Acta Astro 2012; 71: 83–91.

21.

Wei

Lin

Nie

. Investigation on soft-landing dynamics of four-legged lunar lander. Acta Astro 2014; 101: 55–66.

22.

Doiron

Zupp

Jr . Apollo lunar module landing dynamics. In: 41st Structures, Structural Dynamics, and Materials Conference and Exhibit, Atlanta, GA, United States, 3–6 April 2000. AIAA. DOI: 10.2514/6.2000-1678.

23.

Machado

Moreira

Flores

. Compliant contact force models in multibody dynamics: evolution of the Hertz contact theory. Mech Mach Theory 2012; 53: 99–121.

24.

Yang

Zeng

Man

. Design and verification of the landing impact attenuation system for Chang’E-3 lander. Sci Sin (Technologica) [in Chinese] 2014; 44(5): 440–449.

25.

Fang

Sudjianto

. Design and modeling for computer experiments. London, UK: Chapman & Hall/CRC, 2006.

26.

Morris

Mitchell

. Exploratory designs for computational experiments. J Stat Plan Inference 1995; 43(3): 381–402.

27.

Shayanfar

Barkhordari

MA,

Roudak

. An efficient reliability algorithm for locating design point using the combination of importance sampling concepts and response surface method. Commun Nolinear Sci Numer Simul 2017; 47: 223–237.

28.

Draper

Smith

. Applied regression analysis. 3rd ed. New York, United States: John Wiley & Sons, 1998.