#### 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.

#### 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 Figure 1. Since the decoupled equations can be solved by computer directly, the main advantage of modal method is that its computational efficiency is very high.

##### 2.1. Modal Superposition

Continuum body can be described by its mass matrices and stiffness matrices , with a damping matrix ; then the general vibration differential equation is as where represents the displacements vector of the node degrees of freedom (DOF), represents the force vector, and generally has a large amount of DOF . The modal coordinate transformation as shown in (2) is used to decouple vibration differential equations into linear equations: where 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 : in which where is the natural frequency of the th mode , is identity matrix, and and 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).

##### 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: where is the state vector as (6) and is the derivative of the state vector. , , , and , respectively, are the system matrix, input matrix, output matrix, and transfer matrix; and are the I/O of system:

Then (3) is transformed as follows: where is defined in (8); 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: where matrix is used to define the system output. Then each matrix of state space can be determined according (7):

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), Solving the transfer functions 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 on the flexible body is as follows: where is the local position vector, is the transformation matrix, is the local position of point in the undeformed state, and is the deformation vector of point which can be approximated in the modal coordinate using (2).

Thus, we can define the coordinates of body as where and are the reference coordinates respect displacement and rotation, respectively; 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 under the constraints where represents mass submatrices, is the stiffness matrix, is the vector of Lagrange multiplier, matrix represents the transposition of the constraint Jacobian matrix, is the vector of external force, and 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: where is the kinetic energy of body ; and are, respectively, the mass density and volume of body ; and is the global velocity vector of point on the body. The expression of is where is used, the central term on the right-hand side of (18) can in general be written as where is the vector of the time derivatives of the rotational coordinates of the body reference and is defined as where is the total number of rotation coordinates of the reference of body . So (18) can be written as

Equation (17) can be written in a more compact form as

The expression of mass matrix of body is

The mass matrix can also be written as

As given by (23) and (24), the submatrix associated with the translations of the origin of the body reference is a constant matrix; the submatrices and represent the coupling between reference motion and elastic deformation; and submatrices and depend on both the rotational reference coordinates and the elastic coordinates of the body . As (20) shows

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 matrix is the null matrix which means that the matrix is also the null matrix. This, however, is not the case when deformable bodies are considered. As (25) shows, there is no guarantee that is the null matrix because of the time-dependent coordinates . Then vector must be iteratively updated. Therefore, submatrices , , , and must be iteratively updated, respectively. Finally, the submatrix 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:

#### 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 Figure 4 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.

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 Figure 5, where - means the response of vertical when vertical input and - means the response of lateral when vertical input. The magnitude of lateral response ( direction) is greater than that of vertical response ( 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 -.

The 7th mode has the largest frequency magnitude, and, as shown in Figure 6, 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 10^{3} 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.

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 Figure 7(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 where and are, respectively, the vertical stiffness and damping of tire. is the penetration depth.

**(a) Vertical force**

**(b) Lateral force**

The lateral force exerted by the tires is described as a function of slip angle. As shown in Figure 7(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 where is the lateral stiffness factors and is the slip angle. The measured data of tire parameters are shown in Table 3.

The drop multibody system includes leaf spring, load masses, and tires as in Figure 8. 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.

As introduced in Section 3.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 Figure 9. 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.

##### 4.4. Drop Test

As shown in Figure 10(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.

**(a) Drop test system**

**(b) Measuring system**

According to the use-velocity and limit-velocity of landing, set the drop heights as Table 5; the tires are still at the beginning of the drop test.

Video 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 Figure 11. 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.

**(a) Displacement of joint b**

**(b) Reaction force of joint b**

**(c) Displacement of joint d**

**(d) Reaction force of joint d**

In 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 Table 6. 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.

As shown in Table 6 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.

#### 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.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The project was supported by the Knowledge Innovation Program of the Chinese Academy of Sciences (Grant no. YYYJ-1122) and the Innovation Program of UAV funded by the Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Sciences,(CIOMP). The authors also gratefully acknowledge the support from the Science Fund for Young Scholars of the National Natural Science Foundation of China (Grant no. 51305421).

#### Supplementary Materials

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 nonlinear 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.