# A Simplified Flexible Multibody Dynamics for a Main Landing Gear with Flexible Leaf Spring

Authors: Zhi-Peng Xue; Ming Li; Yan-Hui Li; Hong-Guang Jia
Journal: Shock and Vibration
 (2014)

---

## Abstract

The dynamics of multibody systems with deformable components has been a subject of interest in many different fields such as machine design and aerospace. Traditional rigid-flexible systems often take a lot of computer resources to get accurate results. Accuracy and efficiency of computation have been the focus of this research in satisfying the coupling of rigid body and flex body. The method is based on modal analysis and linear theory of elastodynamics: reduced modal datum was used to describe the elastic deformation which was a linear approximate of the flexible part. Then rigid-flexible multibody system was built and the highly nonlinearity of the mass matrix caused by the limited rotation of the deformation part was approximated using the linear theory of elastodynamics. The above methods were used to establish the drop system of the leaf spring type landing gear of a small UAV. Comparisons of the drop test and simulation were applied. Results show that the errors caused by the linear approximation are
acceptable, and the simulation process is fast and stable.


---

## Body

## 1. Introduction

Automatic takeoff, landing, and taxiing are very important parts of completely autonomous flight of UAV, and taxiing tests on UAVs are obviously expensive and risky as it is easy to bring about accidents under high-speed taxiing. Simulation on the other hand is particularly important before the prototype taxiing experiment. During the development of control systems, real-time simulations such as man-in-the-loop (MIL) and hardware-in-the-loop (HIL) are also used to take place between design level simulations and costly experiments with the real plant [1–3]. Obviously accurate and efficient computer simulations of UAV on-ground models play an important role in the control system design, performance evaluations, and dynamics analysis of such vehicle systems. Meanwhile research on ground dynamics in the past decades has achieved significant results, and the design of these dynamics traditionally relies on mathematical model. Ro [4] has studied aircraft-runway dynamics in detail. In general mathematical equations of ground dynamics are complex, highly coupled, and nonlinear; their derivation and numerical implementation demand considerable time and computer cost. In recent years, the accuracy and computational efficiency of flexible multibody has been the focus of research [5, 6]. The simplified flexible multibody method enables the dynamic response to simulate effectively, even be used in many real-time simulations [7–9].Nonretractable type landing gear has been widely used in small UAVs, a leaf spring is used to provide stiffness and damping which absorbs shock and energy by its space deformations [10, 11]. A common way to model body flexibility is the finite element (FE) method which is based on body geometry and can solve deformations stress and strain precisely [12, 13]. But the finite element method often produces a large number of degrees of freedom (DOF) which leads to the computer computation very costly [14]. To reduce the amount of DOF a simplified method was applied which is called modal method that a body’s deformation can be approximated by a linear combination of preselected mode shape functions [15]. Hauser et al. [7] have studied the computational efficiency and error of modal method, and it has been applied in video games, surgery training, and other simulation systems. In this paper modal superposition method was used to establish the flexible model of the leaf spring landing gear.Flexible body, experiencing large space motion while undergoing small elastic deformations, theoretically will cause the inertial coupling between elastic deformation and rigid body motion; the most common approach is to use a time-varying mass matrix which is highly nonlinear and needs to be recomputed in every time step. Many methods [16] are used to simplify the mass matrix. A more simple way is called the linear theory of elastodynamics that the impact of elastic deformation at the multibody dynamics is ignored and may not be suitable for some cases like high-speed rotation.Tire is another elastic part of the system. The most representative models are based on empirical functions, designed and tuned to fit measure data. Some mathematical formulations can be already found in recent research. Pacejka’s Magic Formula (trigonometric formulation) [17] is a popular model in the automotive industry. The NASA, through an AGARD project, proposed a mathematical formulation applied to commercial airplanes and based on polynomial equations [3, 18]. Linear parameter varying (LPV) and linear fractional transformation (LFT) modeling are also used together with a novel force identification approach based on nonlinear dynamic inversion theory to obtain a simplified tire model [19]. To decrease the complexity of the model a typical, widely used tire model is a point contact spring model with either linear or nonlinear stiffness and damping in parallel [4].The idea which motivates this work is to propose a general modeling method for elastic leaf spring. This paper also proposes a simplified flexible multibody to study the dynamics of a leaf spring landing gear, aiming at efficient formulation and solution of equations where a compromise between accuracy and computational efficiency has to be reached. Another contribution of this paper is to discuss the computer implementation of the proposed formulation and demonstrate its use by developing a drop system model. Drop test and simulation are conducted to analyze the model accuracy and efficiency.


## 2. Flexible Model

Modal method is a linear approximation of elastic deformation; complex nonlinear differential equations are decoupled into a set of single DOF linear equations by modal coordinate transformation. These equations can be further transformed into a more briefly state space form as shown in Figure1. Since the decoupled equations can be solved by computer directly, the main advantage of modal method is that its computational efficiency is very high.Figure 1
Block diagram of modal analysis method.


### 2.1. Modal Superposition

Continuum body can be described by its mass matricesM
u and stiffness matrices K
u, with a damping matrix C
u; then the general vibration differential equation is as
(1)
M
u
u
¨
f
+
C
u
u
˙
f
+
K
u
u
f
=
f
,
where u
f represents the displacements vector of the node degrees of freedom (DOF), f represents the force vector, and generally u
f has a large amount of DOF n. The modal coordinate transformation as shown in (2) is used to decouple vibration differential equations into linear equations:
(2)
u
f
=
Φ
q
f
,
where q
f represents the modal coordinates and Ф represents the reduced modal matrix. By substituting (2) into (1) and premultiplying the transpose of matrix Ф, one gets the following reduced system with the dimension n
m
(
n
m
<
n
):
(3)
M
q
q
¨
f
+
C
q
q
˙
f
+
K
q
q
f
=
Q
in which
(4)
M
q
=
Φ
T
M
Φ
=
I
n
m
×
n
m
,
K
q
=
Φ
T
K
Φ
=
diag
⁡
(
ω
i
2
)
,
C
q
=
Φ
T
C
Φ
,
Q
=
Φ
T
f
,
where ω
i is the natural frequency of the ith mode (
1
≤
i
th
≤
n
m
), M
q is identity matrix, and K
q and C
q are diagonal matrices. Thus, the vibration differential equations are approximated by (3). Each modal coordinate can be calculated independently; then modal coordinates (Figure 2) are used to solve for node displacements using (2).Figure 2
Coordinates of flexible body.


### 2.2. State-Space Form

The mathematical calculations through matrix-vector form by computer will be very efficient by transforming (3) into the state space form. The linear state space is as follows:
(5)
x
˙
=
A
x
+
B
u
,
y
=
C
x
+
D
u
,
where x is the state vector as (6) and x
˙ is the derivative of the state vector. A, B, C, and D, respectively, are the system matrix, input matrix, output matrix, and transfer matrix; u and y are the I/O of system:
(6)
x
=
[
q
f
q
˙
f
]
.Then (3) is transformed as follows:
(7)
q
¨
f
+
2
Γ
q
˙
f
+
Ω
q
f
=
Σ
u
,
(8)
Σ
=
Φ
T
F
,
where Σ is defined in (8); F is used to define the degrees of freedom of the system input.According to (5), the displacements, velocities, and accelerations are described by the modal coordinates:
(9)
y
=
[
Θ
q
f
Θ
q
˙
f
Θ
q
¨
f
]
,
Θ
=
T
Φ
,
where matrix T is used to define the system output. Then each matrix of state space can be determined according (7):
(10)
A
=
[
0
I
-
Ω
-
2
Γ
]
,
B
=
[
0
Σ
]
,
C
=
[
Θ
0
0
Θ
-
Θ
Ω
-
2
Θ
Γ
]
,
D
=
[
0
0
Θ
Σ
]
.The matrices of state space are constructed by the reduced modal data, its output is the deformation, and input is the stress of the flexible body.


### 2.3. Modal Reduction

One gets the same number of modes with the number of DOF by modal analysis, and the accuracy of modal method depends on a user-defined frequency range, though the more modal modes are contained, the greater cost of calculation is required. Modal reduction is based on the frequency response analysis which is used to determine the contribution of each mode to the displacement response.Taking Laplace transforms with zero initial conditions of (5),
(11)
x
(
s
)
=
(
s
I
-
A
)
-
1
B
u
(
s
)
,
y
(
s
)
=
C
(
s
I
-
A
)
-
1
B
u
(
s
)
+
D
u
(
s
)
.
Solving the transfer functions
(12)
y
(
s
)
u
(
s
)
=
C
(
s
I
-
A
)
-
1
B
+
D
for combinations of degrees of freedom where forces are applied and where displacements are taken can get the amplitude of the response of each mode.


## 2.1. Modal Superposition

Continuum body can be described by its mass matricesM
u and stiffness matrices K
u, with a damping matrix C
u; then the general vibration differential equation is as
(1)
M
u
u
¨
f
+
C
u
u
˙
f
+
K
u
u
f
=
f
,
where u
f represents the displacements vector of the node degrees of freedom (DOF), f represents the force vector, and generally u
f has a large amount of DOF n. The modal coordinate transformation as shown in (2) is used to decouple vibration differential equations into linear equations:
(2)
u
f
=
Φ
q
f
,
where q
f represents the modal coordinates and Ф represents the reduced modal matrix. By substituting (2) into (1) and premultiplying the transpose of matrix Ф, one gets the following reduced system with the dimension n
m
(
n
m
<
n
):
(3)
M
q
q
¨
f
+
C
q
q
˙
f
+
K
q
q
f
=
Q
in which
(4)
M
q
=
Φ
T
M
Φ
=
I
n
m
×
n
m
,
K
q
=
Φ
T
K
Φ
=
diag
⁡
(
ω
i
2
)
,
C
q
=
Φ
T
C
Φ
,
Q
=
Φ
T
f
,
where ω
i is the natural frequency of the ith mode (
1
≤
i
th
≤
n
m
), M
q is identity matrix, and K
q and C
q are diagonal matrices. Thus, the vibration differential equations are approximated by (3). Each modal coordinate can be calculated independently; then modal coordinates (Figure 2) are used to solve for node displacements using (2).Figure 2
Coordinates of flexible body.


## 2.2. State-Space Form

The mathematical calculations through matrix-vector form by computer will be very efficient by transforming (3) into the state space form. The linear state space is as follows:
(5)
x
˙
=
A
x
+
B
u
,
y
=
C
x
+
D
u
,
where x is the state vector as (6) and x
˙ is the derivative of the state vector. A, B, C, and D, respectively, are the system matrix, input matrix, output matrix, and transfer matrix; u and y are the I/O of system:
(6)
x
=
[
q
f
q
˙
f
]
.Then (3) is transformed as follows:
(7)
q
¨
f
+
2
Γ
q
˙
f
+
Ω
q
f
=
Σ
u
,
(8)
Σ
=
Φ
T
F
,
where Σ is defined in (8); F is used to define the degrees of freedom of the system input.According to (5), the displacements, velocities, and accelerations are described by the modal coordinates:
(9)
y
=
[
Θ
q
f
Θ
q
˙
f
Θ
q
¨
f
]
,
Θ
=
T
Φ
,
where matrix T is used to define the system output. Then each matrix of state space can be determined according (7):
(10)
A
=
[
0
I
-
Ω
-
2
Γ
]
,
B
=
[
0
Σ
]
,
C
=
[
Θ
0
0
Θ
-
Θ
Ω
-
2
Θ
Γ
]
,
D
=
[
0
0
Θ
Σ
]
.The matrices of state space are constructed by the reduced modal data, its output is the deformation, and input is the stress of the flexible body.


## 2.3. Modal Reduction

One gets the same number of modes with the number of DOF by modal analysis, and the accuracy of modal method depends on a user-defined frequency range, though the more modal modes are contained, the greater cost of calculation is required. Modal reduction is based on the frequency response analysis which is used to determine the contribution of each mode to the displacement response.Taking Laplace transforms with zero initial conditions of (5),
(11)
x
(
s
)
=
(
s
I
-
A
)
-
1
B
u
(
s
)
,
y
(
s
)
=
C
(
s
I
-
A
)
-
1
B
u
(
s
)
+
D
u
(
s
)
.
Solving the transfer functions
(12)
y
(
s
)
u
(
s
)
=
C
(
s
I
-
A
)
-
1
B
+
D
for combinations of degrees of freedom where forces are applied and where displacements are taken can get the amplitude of the response of each mode.


## 3. Rigid-Flexible Multibody Dynamics

### 3.1. Multibody Dynamics

The most frequent approach used in multibody system to consider flexibility is the floating frame of reference formulation (FFRF), in which two sets of coordinate are used to describe a deformable body as discussed by Shabana [20]. One set is reference coordinate which describes the location and orientation of a selected body coordinate system, while the second one is elastic coordinate which describes the deformation of the body with respect to its coordinates. The global position of an arbitrary point B on the flexible body is as follows:
(13)
r
i
B
=
r
i
+
A
i
(
u
i
B
+
u
f
B
)
,
where r
i is the local position vector, A
i is the transformation matrix, u
i
B is the local position of point B in the undeformed state, and u
f
B is the deformation vector of point B which can be approximated in the modal coordinate using (2).Thus, we can define the coordinates of bodyi as
(14)
q
i
=
[
R
i
θ
i
q
f
i
]
,
where R
i and θ
i are the reference coordinates respect displacement and rotation, respectively; q
f
i is the vector of elastic coordinates.The equations of motion of Lagrangian formalism in multibody systems have been proposed in [20]. The augmented matrix equation expressed in generalized coordinates is
(15)
[
m
R
R
i
m
R
θ
i
m
R
f
i
m
θ
θ
i
m
θ
f
i
symmetric
m
f
f
i
]
[
R
¨
i
θ
¨
i
q
¨
f
i
]
+
[
0
0
0
0
0
0
0
0
K
f
f
i
]
[
R
i
θ
i
q
f
i
]
+
[
C
R
i
T
C
θ
i
T
C
q
f
i
T
]
λ
=
[
(
Q
e
i
)
R
(
Q
e
i
)
θ
(
Q
e
i
)
f
]
+
[
(
Q
v
i
)
R
(
Q
v
i
)
θ
(
Q
v
i
)
f
]
under the constraints
(16)
C
(
R
,
θ
,
q
f
)
=
0
,
where m
i
i represents mass submatrices,  K
f
f is the stiffness matrix, λ is the vector of Lagrange multiplier, matrix C
q
T represents the transposition of the constraint Jacobian matrix, Q
e is the vector of external force, and Q
v is the vector of quadratic velocity.


### 3.2. Coupling between Reference and Elastic Displacements

#### 3.2.1. Inertia of Deformable Bodies

The following definition of the kinetic energy is used:(17)
T
i
=
1
2
∫
V
i
ρ
i
r
˙
i
B
T
r
˙
i
B
d
V
i
,
where T
i is the kinetic energy of body i; ρ
i and V
i are, respectively, the mass density and volume of body i; and r
˙
i
B is the global velocity vector of point B on the body. The expression of r
˙
i
B is
(18)
r
˙
i
B
=
r
˙
i
+
A
˙
i
(
u
i
B
+
u
f
B
)
+
A
i
(
u
˙
i
B
+
u
˙
f
B
)
,
where u
˙
i
B
=
0 is used, the central term on the right-hand side of (18) can in general be written as
(19)
A
˙
i
(
u
i
B
+
u
f
B
)
=
B
i
θ
˙
i
,
where θ
˙
i is the vector of the time derivatives of the rotational coordinates of the body reference and B
i is defined as
(20)
B
i
=
[
∂
∂
θ
1
i
A
i
(
u
i
B
+
u
f
B
)
⋯
∂
∂
θ
n
r
i
A
i
(
u
i
B
+
u
f
B
)
]
,
where n
r is the total number of rotation coordinates of the reference of body i. So (18) can be written as
(21)
r
˙
i
B
=
[
I
B
i
A
i
Φ
i
]
[
r
˙
i
θ
˙
i
q
˙
f
i
]
.Equation (17) can be written in a more compact form as
(22)
T
i
=
1
2
q
˙
T
i
M
i
q
˙
i
.The expression of mass matrix of bodyi is
(23)
M
i
=
∫
V
i
ρ
i
[
I
B
i
T
(
A
i
Φ
i
)
T
]
[
I
B
i
A
i
Φ
i
]
d
V
i
=
∫
V
i
ρ
i
[
I
B
i
A
i
Φ
i
B
i
T
B
i
B
i
T
A
i
Φ
i
symmetric
Φ
i
T
Φ
i
]
d
V
i
.The mass matrix can also be written as(24)
M
i
=
[
m
R
R
i
m
R
θ
i
m
R
f
i
m
θ
θ
i
m
θ
f
i
symmetric
m
f
f
i
]
.As given by (23) and (24), the submatrix m
R
R
i associated with the translations of the origin of the body reference is a constant matrix; the submatrices m
R
f
i and m
θ
f
i represent the coupling between reference motion and elastic deformation; and submatrices m
θ
θ
i and m
R
θ
i depend on both the rotational reference coordinates and the elastic coordinates of the body i. As (20) shows
(25)
B
i
=
B
i
(
θ
˙
i
,
q
˙
f
i
)
.If the body is rigid and the origin of the body reference is attached to the center of mass, the translation and rotation of the rigid body are dynamically decoupled. In this case the matrixm
R
θ
i is the null matrix which means that the matrix B
i is also the null matrix. This, however, is not the case when deformable bodies are considered. As (25) shows, there is no guarantee that B
i is the null matrix because of the time-dependent coordinates q
˙
f
i. Then vector B
i must be iteratively updated. Therefore, submatrices m
R
θ
i, m
θ
f
i, m
R
f
i, and m
θ
θ
i must be iteratively updated, respectively. Finally, the submatrix m
f
f
i in (24) is independent of the generalized coordinates of the body, therefore, is constant.


#### 3.2.2. Linear Theory of Elastodynamics

As above section shows, the mixed set of reference and elastic coordinates of deformable body leads to a highly nonlinear mass matrix as a result of the inertia coupling between the reference and the elastic displacements. The mass matrix must be iteratively updated.Rigid body motion and flexible deformation interact with each other though inertial effects. The deformation body is rotating rapidly; neither effect will be significant. Therefore a solution strategy that has been used to omit these effects. The multibody system was first as rigid body systems as (26). General multi-rigid-body computer programs then be used to solve for the inertia and reaction forces. These forces are then introduced to the state space form of elastic deformable model to solve the deformation displacements. The total motion of a body is then obtained by superimposing the small elastic deformation on the gross rigid body motion as shown in Figure 3:
(26)
[
m
R
R
i
m
θ
θ
i
]
[
R
¨
i
θ
¨
i
]
+
[
C
R
i
T
C
θ
i
T
]
λ
=
[
(
Q
e
i
)
R
(
Q
e
i
)
θ
]
+
[
(
Q
v
i
)
R
(
Q
v
i
)
θ
]
.Figure 3
Block diagram of rigid-flexible multibody system.


## 3.1. Multibody Dynamics

The most frequent approach used in multibody system to consider flexibility is the floating frame of reference formulation (FFRF), in which two sets of coordinate are used to describe a deformable body as discussed by Shabana [20]. One set is reference coordinate which describes the location and orientation of a selected body coordinate system, while the second one is elastic coordinate which describes the deformation of the body with respect to its coordinates. The global position of an arbitrary point B on the flexible body is as follows:
(13)
r
i
B
=
r
i
+
A
i
(
u
i
B
+
u
f
B
)
,
where r
i is the local position vector, A
i is the transformation matrix, u
i
B is the local position of point B in the undeformed state, and u
f
B is the deformation vector of point B which can be approximated in the modal coordinate using (2).Thus, we can define the coordinates of bodyi as
(14)
q
i
=
[
R
i
θ
i
q
f
i
]
,
where R
i and θ
i are the reference coordinates respect displacement and rotation, respectively; q
f
i is the vector of elastic coordinates.The equations of motion of Lagrangian formalism in multibody systems have been proposed in [20]. The augmented matrix equation expressed in generalized coordinates is
(15)
[
m
R
R
i
m
R
θ
i
m
R
f
i
m
θ
θ
i
m
θ
f
i
symmetric
m
f
f
i
]
[
R
¨
i
θ
¨
i
q
¨
f
i
]
+
[
0
0
0
0
0
0
0
0
K
f
f
i
]
[
R
i
θ
i
q
f
i
]
+
[
C
R
i
T
C
θ
i
T
C
q
f
i
T
]
λ
=
[
(
Q
e
i
)
R
(
Q
e
i
)
θ
(
Q
e
i
)
f
]
+
[
(
Q
v
i
)
R
(
Q
v
i
)
θ
(
Q
v
i
)
f
]
under the constraints
(16)
C
(
R
,
θ
,
q
f
)
=
0
,
where m
i
i represents mass submatrices,  K
f
f is the stiffness matrix, λ is the vector of Lagrange multiplier, matrix C
q
T represents the transposition of the constraint Jacobian matrix, Q
e is the vector of external force, and Q
v is the vector of quadratic velocity.


## 3.2. Coupling between Reference and Elastic Displacements

### 3.2.1. Inertia of Deformable Bodies

The following definition of the kinetic energy is used:(17)
T
i
=
1
2
∫
V
i
ρ
i
r
˙
i
B
T
r
˙
i
B
d
V
i
,
where T
i is the kinetic energy of body i; ρ
i and V
i are, respectively, the mass density and volume of body i; and r
˙
i
B is the global velocity vector of point B on the body. The expression of r
˙
i
B is
(18)
r
˙
i
B
=
r
˙
i
+
A
˙
i
(
u
i
B
+
u
f
B
)
+
A
i
(
u
˙
i
B
+
u
˙
f
B
)
,
where u
˙
i
B
=
0 is used, the central term on the right-hand side of (18) can in general be written as
(19)
A
˙
i
(
u
i
B
+
u
f
B
)
=
B
i
θ
˙
i
,
where θ
˙
i is the vector of the time derivatives of the rotational coordinates of the body reference and B
i is defined as
(20)
B
i
=
[
∂
∂
θ
1
i
A
i
(
u
i
B
+
u
f
B
)
⋯
∂
∂
θ
n
r
i
A
i
(
u
i
B
+
u
f
B
)
]
,
where n
r is the total number of rotation coordinates of the reference of body i. So (18) can be written as
(21)
r
˙
i
B
=
[
I
B
i
A
i
Φ
i
]
[
r
˙
i
θ
˙
i
q
˙
f
i
]
.Equation (17) can be written in a more compact form as
(22)
T
i
=
1
2
q
˙
T
i
M
i
q
˙
i
.The expression of mass matrix of bodyi is
(23)
M
i
=
∫
V
i
ρ
i
[
I
B
i
T
(
A
i
Φ
i
)
T
]
[
I
B
i
A
i
Φ
i
]
d
V
i
=
∫
V
i
ρ
i
[
I
B
i
A
i
Φ
i
B
i
T
B
i
B
i
T
A
i
Φ
i
symmetric
Φ
i
T
Φ
i
]
d
V
i
.The mass matrix can also be written as(24)
M
i
=
[
m
R
R
i
m
R
θ
i
m
R
f
i
m
θ
θ
i
m
θ
f
i
symmetric
m
f
f
i
]
.As given by (23) and (24), the submatrix m
R
R
i associated with the translations of the origin of the body reference is a constant matrix; the submatrices m
R
f
i and m
θ
f
i represent the coupling between reference motion and elastic deformation; and submatrices m
θ
θ
i and m
R
θ
i depend on both the rotational reference coordinates and the elastic coordinates of the body i. As (20) shows
(25)
B
i
=
B
i
(
θ
˙
i
,
q
˙
f
i
)
.If the body is rigid and the origin of the body reference is attached to the center of mass, the translation and rotation of the rigid body are dynamically decoupled. In this case the matrixm
R
θ
i is the null matrix which means that the matrix B
i is also the null matrix. This, however, is not the case when deformable bodies are considered. As (25) shows, there is no guarantee that B
i is the null matrix because of the time-dependent coordinates q
˙
f
i. Then vector B
i must be iteratively updated. Therefore, submatrices m
R
θ
i, m
θ
f
i, m
R
f
i, and m
θ
θ
i must be iteratively updated, respectively. Finally, the submatrix m
f
f
i in (24) is independent of the generalized coordinates of the body, therefore, is constant.


### 3.2.2. Linear Theory of Elastodynamics

As above section shows, the mixed set of reference and elastic coordinates of deformable body leads to a highly nonlinear mass matrix as a result of the inertia coupling between the reference and the elastic displacements. The mass matrix must be iteratively updated.Rigid body motion and flexible deformation interact with each other though inertial effects. The deformation body is rotating rapidly; neither effect will be significant. Therefore a solution strategy that has been used to omit these effects. The multibody system was first as rigid body systems as (26). General multi-rigid-body computer programs then be used to solve for the inertia and reaction forces. These forces are then introduced to the state space form of elastic deformable model to solve the deformation displacements. The total motion of a body is then obtained by superimposing the small elastic deformation on the gross rigid body motion as shown in Figure 3:
(26)
[
m
R
R
i
m
θ
θ
i
]
[
R
¨
i
θ
¨
i
]
+
[
C
R
i
T
C
θ
i
T
]
λ
=
[
(
Q
e
i
)
R
(
Q
e
i
)
θ
]
+
[
(
Q
v
i
)
R
(
Q
v
i
)
θ
]
.Figure 3
Block diagram of rigid-flexible multibody system.


## 3.2.1. Inertia of Deformable Bodies

The following definition of the kinetic energy is used:(17)
T
i
=
1
2
∫
V
i
ρ
i
r
˙
i
B
T
r
˙
i
B
d
V
i
,
where T
i is the kinetic energy of body i; ρ
i and V
i are, respectively, the mass density and volume of body i; and r
˙
i
B is the global velocity vector of point B on the body. The expression of r
˙
i
B is
(18)
r
˙
i
B
=
r
˙
i
+
A
˙
i
(
u
i
B
+
u
f
B
)
+
A
i
(
u
˙
i
B
+
u
˙
f
B
)
,
where u
˙
i
B
=
0 is used, the central term on the right-hand side of (18) can in general be written as
(19)
A
˙
i
(
u
i
B
+
u
f
B
)
=
B
i
θ
˙
i
,
where θ
˙
i is the vector of the time derivatives of the rotational coordinates of the body reference and B
i is defined as
(20)
B
i
=
[
∂
∂
θ
1
i
A
i
(
u
i
B
+
u
f
B
)
⋯
∂
∂
θ
n
r
i
A
i
(
u
i
B
+
u
f
B
)
]
,
where n
r is the total number of rotation coordinates of the reference of body i. So (18) can be written as
(21)
r
˙
i
B
=
[
I
B
i
A
i
Φ
i
]
[
r
˙
i
θ
˙
i
q
˙
f
i
]
.Equation (17) can be written in a more compact form as
(22)
T
i
=
1
2
q
˙
T
i
M
i
q
˙
i
.The expression of mass matrix of bodyi is
(23)
M
i
=
∫
V
i
ρ
i
[
I
B
i
T
(
A
i
Φ
i
)
T
]
[
I
B
i
A
i
Φ
i
]
d
V
i
=
∫
V
i
ρ
i
[
I
B
i
A
i
Φ
i
B
i
T
B
i
B
i
T
A
i
Φ
i
symmetric
Φ
i
T
Φ
i
]
d
V
i
.The mass matrix can also be written as(24)
M
i
=
[
m
R
R
i
m
R
θ
i
m
R
f
i
m
θ
θ
i
m
θ
f
i
symmetric
m
f
f
i
]
.As given by (23) and (24), the submatrix m
R
R
i associated with the translations of the origin of the body reference is a constant matrix; the submatrices m
R
f
i and m
θ
f
i represent the coupling between reference motion and elastic deformation; and submatrices m
θ
θ
i and m
R
θ
i depend on both the rotational reference coordinates and the elastic coordinates of the body i. As (20) shows
(25)
B
i
=
B
i
(
θ
˙
i
,
q
˙
f
i
)
.If the body is rigid and the origin of the body reference is attached to the center of mass, the translation and rotation of the rigid body are dynamically decoupled. In this case the matrixm
R
θ
i is the null matrix which means that the matrix B
i is also the null matrix. This, however, is not the case when deformable bodies are considered. As (25) shows, there is no guarantee that B
i is the null matrix because of the time-dependent coordinates q
˙
f
i. Then vector B
i must be iteratively updated. Therefore, submatrices m
R
θ
i, m
θ
f
i, m
R
f
i, and m
θ
θ
i must be iteratively updated, respectively. Finally, the submatrix m
f
f
i in (24) is independent of the generalized coordinates of the body, therefore, is constant.


## 3.2.2. Linear Theory of Elastodynamics

As above section shows, the mixed set of reference and elastic coordinates of deformable body leads to a highly nonlinear mass matrix as a result of the inertia coupling between the reference and the elastic displacements. The mass matrix must be iteratively updated.Rigid body motion and flexible deformation interact with each other though inertial effects. The deformation body is rotating rapidly; neither effect will be significant. Therefore a solution strategy that has been used to omit these effects. The multibody system was first as rigid body systems as (26). General multi-rigid-body computer programs then be used to solve for the inertia and reaction forces. These forces are then introduced to the state space form of elastic deformable model to solve the deformation displacements. The total motion of a body is then obtained by superimposing the small elastic deformation on the gross rigid body motion as shown in Figure 3:
(26)
[
m
R
R
i
m
θ
θ
i
]
[
R
¨
i
θ
¨
i
]
+
[
C
R
i
T
C
θ
i
T
]
λ
=
[
(
Q
e
i
)
R
(
Q
e
i
)
θ
]
+
[
(
Q
v
i
)
R
(
Q
v
i
)
θ
]
.Figure 3
Block diagram of rigid-flexible multibody system.


## 4. Simulation and Test

Sections above have presented a very effective modeling method for rigid-flexible multibody system which will be used in this section to build the main landing gear of a small UAV in Matlab/Simulink. Tests are implemented to verify the accuracy of this method.


### 4.1. Modeling of the Leaf Spring

Cantilever type leaf spring was manufactured by composite material layers as shown in Figure4 and the parameters of each material are shown in Table 1. The FE model of leaf spring has 2868 nodes which means that it has 2868 × 3 modes. The first 6 modes are rigid body modes resulting from a free-free modal analysis and the 7th to 17th modes are shown in Table 2.Table 1
Layer parameters.
Parameters
Materials
Carbon cloth
Glass fibre
Glass prepreg
Thickness (mm)
0.2
0.2
0.2
Young’s modulusE
x (Gpa)
57–66
17–20
42–52
Young’s modulusE
y (Gpa)
57–66
17–20
8–10
Young’s modulusE
z (Gpa)
8–10
8–10
8–10
Shear modulusG (Gpa)
4.8
4
3.5
Poisson’s parameters
0.3
0.3
0.3
Destiny (g/cm3)
1.6
2
2Table 2
Modal analysis results.
Modes
Frequency rad/s
7
107
8
322
9
494
10
641
11
1261
12
1294
13
1927
14
2945
15
3393
16
3453
17
3849Figure 4
Five-layer leaf spring landing gear.Frequency response analysis for first 50 modes except the rigid body modes is performed at the joint connected the landing gear and tire. The results are shown in Figure5, where Z-Z means the response of vertical when vertical input and Z-Y means the response of lateral when vertical input. The magnitude of lateral response (Y direction) is greater than that of vertical response (Z direction) at the former modes, which means that the lateral stiffness is lower than the vertical stiffness. For the reduced modes must include the major modes of the two directions, mode reduce is based on Z-Z.Figure 5
Frequency response.The 7th mode has the largest frequency magnitude, and, as shown in Figure6, 7th mode shape is the first bending mode; this bending shape was appearing repeatedly during the compressing of the leaf spring of drop tests, which indicated that the 7th mode is the most important mode. The magnitude response decreases as the frequency increases, and the first area peak appears at the 10th mode, and its mode shape is second bending mode. The mode which has a resonant area peak is the important mode, such as the 11th mode; it is the third bending mode. The magnitude of the 17th mode is about 103 times smaller than that of the 7th mode, and the magnitude of the resonant peaks after 17th modes is far smaller than that of the 17th mode. So the 7th to 17th modes are the major modes.Figure 6
Parts of mode shapes.Then 11 modes are used to construct the elastic model, and each mode shape of the reduced modes still contains the modal datum of every node. For we only concern the nodes where forces are applied and where displacements are taken, the nodes except these located on the joints and where we want to measure the deformations are excluded. So the flexible model is constructed by the reduced modes and reduced nodes.


### 4.2. Rigid-Flexible Multibody Model

A penetration formulation is used to calculate the normal force between the tire and ground as shown in Figure7(a). Similar to a spring/damper formula, the force applied to the tire is proportional to the penetration depth and the penetration rate. The normal force is described as
(27)
F
nom
⁡
=
K
nom
⁡
τ
+
C
nom
⁡
τ
˙
,
where K
nom
⁡ and C
nom
⁡ are, respectively, the vertical stiffness and damping of tire. τ is the penetration depth.Tire model.
(a)
Vertical force
(b)
Lateral forceThe lateral force exerted by the tires is described as a function of slip angle. As shown in Figure7(b), slip angle β is the orientations of the velocity of the tire’s points of contact the surface of motion with respect to the tires’ longitudinal axes. The lateral force is applied to account for the side loading on the tire. It is applied in a direction aligned with the rotation axis of the wheel. This force acts in the opposite direction to the lateral motion of the wheel and can be written as
(28)
F
lat
=
C
1
β
,
where C
1 is the lateral stiffness factors and β is the slip angle. The measured data of tire parameters are shown in Table 3.Table 3
Tire parameters.
Parameter
Value
Radius
0.1 m
Width
0.075 m
Vertical stiffness
75000 N/m
Lateral stiffness
166.55 N/°
Mass
1.5 kg
I
x
x
/
I
y
y
3.45 * 10−3 Kgm2
I
z
z
5.62 * 10−3 Kgm2
Damping
175 N/(m/s)The drop multibody system includes leaf spring, load masses, and tires as in Figure8. The leaf spring and load masses are connected by 2 pin joints a and b, the leaf spring and tires are connected by 2 revolute joints c and d, and the landing gear system is connected with ground by a prismatic joint.Figure 8
Drop multibody system.As introduced in Section3.2, both rigid model and flexible model are used to build the leaf spring, and the information exchanged between them is reaction forces and deformations; the external forces of the flexible model are the reaction forces of the joints. Massless bodies are used to transfer information between rigid body and flexible body. The entire model experiences free falling if the tires have not contacted the ground. If contact happens, reaction forces are solved by rigid system simulation; then these forces are as the inputs of the state space to calculate the displacements of the joint nodes which are used to correct the global positions of these joint nodes in the next simulation step; rigid-flexible multibody simulation is accomplished as this iterative loop.


### 4.3. Statics Test

Statics test was conducted before the drop test as shown in Figure9. As the results show in Table 4, the errors of vertical displacement are below 3% for the statics conditions. So the present flexible model of leaf spring leads to excellent accuracy under the statics situation.Table 4
Statics results.
Load mass
Tests (mm)
Simulations (mm)
Errors
10 kg
5.125
5.079
−0.9%
20 kg
10.250
10.11
−0.6%
30 kg
15.750
15.754
0
40 kg
21.000
21.4145
2%
50 kg
26.625
27.3015
2.5%Figure 9
Statics test.


### 4.4. Drop Test

As shown in Figure10(a), the drop test system includes supporting subsystem, slide track subsystem, lift subsystem, and data acquisition subsystem. The drop mass (100 kg) is attached to the top of the landing gear to simulate the landing mass. For the test results are used to compare with the simulation results, so only the vertical velocity is considered.Devices of drop test system.
(a)
Drop test system
(b)
Measuring systemAccording to the use-velocity and limit-velocity of landing, set the drop heights as Table5; the tires are still at the beginning of the drop test.Table 5
Parameters of drop tests.
Number
Drop velocity
Load mass
Drop height
Drop attitude
1
1.8 m/s
180 kg
0.165 m
Level
2
3 m/s
180 kg
0.46 m
LevelVideo footage is used to capture the vertical and lateral displacements during drop tests by the vertical and lateral fixed rulers.


### 4.5. Results

The flexible body motion of landing gear is shown in Figure11. Results including the displacements and reaction forces of simulation 2 are shown in Figure 12. The simulation step is 1 ms. It is because of the loss of energy during impact which is caused by the frictions of drop test that the shock time of drop test is shorter than that of the drop simulation.Figure 11
Motions of joint b of landing gear.Simulation results of leaf spring of simulation 2.
(a)
Displacement of joint b
(b)
Reaction force of joint b
(c)
Displacement of joint d
(d)
Reaction force of joint dIn test 2 the vertical displacements of joint b from the measuring results of video footage and simulation results when the landing gear undergoing fully compressed are given in Table6. In which a negative displacement means the landing gear is under fully compress and a positive displacement means the landing gear is on the rebound.Table 6
Displacements of joint b in test 2 (mm).
Drop test
Simulation
Error
−155.4
−141.9
8.7%
122.5
132.7
8.3%
−61.3
−64.7
5.5%
26.6
28.7
7.9%
−10.8
−11.4
5.6%As shown in Table6 the landing gear absorbed majority energy after 3 compressions, and the errors are less than 10%. Elastic deformations of joint d during the three compressions are shown in Table 7.Table 7
Error of elastic deformations of joint d in test 2 (mm).
Compression number
Vertical
Lateral
1
8.0%
7.6%
2
6.4%
6.1%
3
5.3%
4.9%


## 4.1. Modeling of the Leaf Spring

Cantilever type leaf spring was manufactured by composite material layers as shown in Figure4 and the parameters of each material are shown in Table 1. The FE model of leaf spring has 2868 nodes which means that it has 2868 × 3 modes. The first 6 modes are rigid body modes resulting from a free-free modal analysis and the 7th to 17th modes are shown in Table 2.Table 1
Layer parameters.
Parameters
Materials
Carbon cloth
Glass fibre
Glass prepreg
Thickness (mm)
0.2
0.2
0.2
Young’s modulusE
x (Gpa)
57–66
17–20
42–52
Young’s modulusE
y (Gpa)
57–66
17–20
8–10
Young’s modulusE
z (Gpa)
8–10
8–10
8–10
Shear modulusG (Gpa)
4.8
4
3.5
Poisson’s parameters
0.3
0.3
0.3
Destiny (g/cm3)
1.6
2
2Table 2
Modal analysis results.
Modes
Frequency rad/s
7
107
8
322
9
494
10
641
11
1261
12
1294
13
1927
14
2945
15
3393
16
3453
17
3849Figure 4
Five-layer leaf spring landing gear.Frequency response analysis for first 50 modes except the rigid body modes is performed at the joint connected the landing gear and tire. The results are shown in Figure5, where Z-Z means the response of vertical when vertical input and Z-Y means the response of lateral when vertical input. The magnitude of lateral response (Y direction) is greater than that of vertical response (Z direction) at the former modes, which means that the lateral stiffness is lower than the vertical stiffness. For the reduced modes must include the major modes of the two directions, mode reduce is based on Z-Z.Figure 5
Frequency response.The 7th mode has the largest frequency magnitude, and, as shown in Figure6, 7th mode shape is the first bending mode; this bending shape was appearing repeatedly during the compressing of the leaf spring of drop tests, which indicated that the 7th mode is the most important mode. The magnitude response decreases as the frequency increases, and the first area peak appears at the 10th mode, and its mode shape is second bending mode. The mode which has a resonant area peak is the important mode, such as the 11th mode; it is the third bending mode. The magnitude of the 17th mode is about 103 times smaller than that of the 7th mode, and the magnitude of the resonant peaks after 17th modes is far smaller than that of the 17th mode. So the 7th to 17th modes are the major modes.Figure 6
Parts of mode shapes.Then 11 modes are used to construct the elastic model, and each mode shape of the reduced modes still contains the modal datum of every node. For we only concern the nodes where forces are applied and where displacements are taken, the nodes except these located on the joints and where we want to measure the deformations are excluded. So the flexible model is constructed by the reduced modes and reduced nodes.


## 4.2. Rigid-Flexible Multibody Model

A penetration formulation is used to calculate the normal force between the tire and ground as shown in Figure7(a). Similar to a spring/damper formula, the force applied to the tire is proportional to the penetration depth and the penetration rate. The normal force is described as
(27)
F
nom
⁡
=
K
nom
⁡
τ
+
C
nom
⁡
τ
˙
,
where K
nom
⁡ and C
nom
⁡ are, respectively, the vertical stiffness and damping of tire. τ is the penetration depth.Tire model.
(a)
Vertical force
(b)
Lateral forceThe lateral force exerted by the tires is described as a function of slip angle. As shown in Figure7(b), slip angle β is the orientations of the velocity of the tire’s points of contact the surface of motion with respect to the tires’ longitudinal axes. The lateral force is applied to account for the side loading on the tire. It is applied in a direction aligned with the rotation axis of the wheel. This force acts in the opposite direction to the lateral motion of the wheel and can be written as
(28)
F
lat
=
C
1
β
,
where C
1 is the lateral stiffness factors and β is the slip angle. The measured data of tire parameters are shown in Table 3.Table 3
Tire parameters.
Parameter
Value
Radius
0.1 m
Width
0.075 m
Vertical stiffness
75000 N/m
Lateral stiffness
166.55 N/°
Mass
1.5 kg
I
x
x
/
I
y
y
3.45 * 10−3 Kgm2
I
z
z
5.62 * 10−3 Kgm2
Damping
175 N/(m/s)The drop multibody system includes leaf spring, load masses, and tires as in Figure8. The leaf spring and load masses are connected by 2 pin joints a and b, the leaf spring and tires are connected by 2 revolute joints c and d, and the landing gear system is connected with ground by a prismatic joint.Figure 8
Drop multibody system.As introduced in Section3.2, both rigid model and flexible model are used to build the leaf spring, and the information exchanged between them is reaction forces and deformations; the external forces of the flexible model are the reaction forces of the joints. Massless bodies are used to transfer information between rigid body and flexible body. The entire model experiences free falling if the tires have not contacted the ground. If contact happens, reaction forces are solved by rigid system simulation; then these forces are as the inputs of the state space to calculate the displacements of the joint nodes which are used to correct the global positions of these joint nodes in the next simulation step; rigid-flexible multibody simulation is accomplished as this iterative loop.


## 4.3. Statics Test

Statics test was conducted before the drop test as shown in Figure9. As the results show in Table 4, the errors of vertical displacement are below 3% for the statics conditions. So the present flexible model of leaf spring leads to excellent accuracy under the statics situation.Table 4
Statics results.
Load mass
Tests (mm)
Simulations (mm)
Errors
10 kg
5.125
5.079
−0.9%
20 kg
10.250
10.11
−0.6%
30 kg
15.750
15.754
0
40 kg
21.000
21.4145
2%
50 kg
26.625
27.3015
2.5%Figure 9
Statics test.


## 4.4. Drop Test

As shown in Figure10(a), the drop test system includes supporting subsystem, slide track subsystem, lift subsystem, and data acquisition subsystem. The drop mass (100 kg) is attached to the top of the landing gear to simulate the landing mass. For the test results are used to compare with the simulation results, so only the vertical velocity is considered.Devices of drop test system.
(a)
Drop test system
(b)
Measuring systemAccording to the use-velocity and limit-velocity of landing, set the drop heights as Table5; the tires are still at the beginning of the drop test.Table 5
Parameters of drop tests.
Number
Drop velocity
Load mass
Drop height
Drop attitude
1
1.8 m/s
180 kg
0.165 m
Level
2
3 m/s
180 kg
0.46 m
LevelVideo footage is used to capture the vertical and lateral displacements during drop tests by the vertical and lateral fixed rulers.


## 4.5. Results

The flexible body motion of landing gear is shown in Figure11. Results including the displacements and reaction forces of simulation 2 are shown in Figure 12. The simulation step is 1 ms. It is because of the loss of energy during impact which is caused by the frictions of drop test that the shock time of drop test is shorter than that of the drop simulation.Figure 11
Motions of joint b of landing gear.Simulation results of leaf spring of simulation 2.
(a)
Displacement of joint b
(b)
Reaction force of joint b
(c)
Displacement of joint d
(d)
Reaction force of joint dIn test 2 the vertical displacements of joint b from the measuring results of video footage and simulation results when the landing gear undergoing fully compressed are given in Table6. In which a negative displacement means the landing gear is under fully compress and a positive displacement means the landing gear is on the rebound.Table 6
Displacements of joint b in test 2 (mm).
Drop test
Simulation
Error
−155.4
−141.9
8.7%
122.5
132.7
8.3%
−61.3
−64.7
5.5%
26.6
28.7
7.9%
−10.8
−11.4
5.6%As shown in Table6 the landing gear absorbed majority energy after 3 compressions, and the errors are less than 10%. Elastic deformations of joint d during the three compressions are shown in Table 7.Table 7
Error of elastic deformations of joint d in test 2 (mm).
Compression number
Vertical
Lateral
1
8.0%
7.6%
2
6.4%
6.1%
3
5.3%
4.9%


## 5. Conclusions

A simplified rigid-flexible multibody method is presented and used in modeling the main landing gear system of a small UAV. The computational efficiency of the vibration differential equations which has been decoupled by modal coordinate transformation is very high. Modal superposition has greatly reduced the model that only 11 modes are used to build the flexible model of the leaf spring, while the FE model has 17208 DOFs. The highly nonlinear mass matrix because of the finite rotation of the flexible part is eliminated by linear elastic dynamics. So the whole model is characterized by a very small number of equations, and there are no significant nonlinear factors that a very high computational efficiency can be reached and it is very suitable for real-time simulation systems.The comparisons between test and simulations are applied to verify the accuracy of this approach. Static errors are less than 5%, and dynamics errors are less than 10%; meanwhile the simulation step has reached to 1 ms. So the computation of this approach is accuracy and efficiency. The modal datum used by this approach can be get from not only the FE modal analysis, but also the modal experiments to modify the first few modes, to get a more accurate and reliable result. Results have shown that the proposed method has a good potential to model UAV with elastic leaf spring landing gear on-ground multibody system which can be used in real-time simulations such as MIL and HIL simulations in the future.



---

*Source: 595964-2014-07-23.xml*