# Dynamic analysis of semi-flexible multibody systems: Coordinates relatives method

Alexis Mouhingou, Naoufel Azouz

<figure><img src="image_1.png" /></figure>

To cite this version:

Alexis Mouhingou, Naoufel Azouz. Dynamic analysis of semi-flexible multibody systems: Coordinates relatives method. ASME International Mechanical Engineering Congress and Exposition (IMECE 2008), Oct 2009, Boston, United States. pp.429-438, 10.1115/IMECE2008-69087. hal-00869654

HAL Id: hal-00869654

https://hal.science/hal-00869654

Submitted on 6 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# DYNAMIC ANALYSIS OF SEMI-FLEXIBLE MULTIBODY SYSTEMS:
## COORDINATES RELATIVES METHOD

MOUHINGOU Alexis, AZOUZ Naoufel

Laboratoire Informatique, Biologie Intégrative et Systèmes Complexes

Université d'Evry Val d'Essonne

40, rue du Pelvoux, 91020 Evry cedex, France

E-mail: Alexis.Mouhingou, Naoufel.Azouz@iup.univ-evry.fr

**Abstract:**

This work presents the dynamic modelling of a multibody systems in cross form constituted of a central body which is connected four flexible arms, at each end of arm is connected a rotor. A particular attention is given to the influence of flexibility on the dynamic behaviour of system. For elaborate the complete dynamic model, one consider the sub-structuration technique in using the Lagrangian approach based on the relatives coordinates method of central body. One establishes then the mathematics equations permitting to obtain the dynamic model of motion for the simulation and the control study. The aerodynamic loads and the gravity force are taking into account for the dynamic complete model. At the end, one considers the case of stationary flight of a miniature Quadrotor. The numerical results permit to simulate the motion of Quadrirotor in flight. But also to compare the flexible and rigid model in order to resort the flexibility effects.

**Keywords:** modelling, flexible multibody systems, Lagrangian dynamics, sub-structuration, UAV, aerodynamic.

**INTRODUCTION:**

The dynamic modelling of a mechanical system in flight presents some difficulties due to the topology of the system but also the take into account of aerodynamic loads, above all when this one is considered as multibody system [3], [4], [12], [13], [14], [17] and [19]. For this purpose, several techniques are proposed in the literature. However the presence of the flexible substructures constitutes an inherent problem of the behaviour study of those machines. For that it is necessary to add the effects of aerodynamic loads in the case of flying objects. It should be noted that, the presence of the various joins between different substructures of the multibody systems make more complex the modelling of such systems. A complete dynamic

model permitting to take into account those different aspects would be a considerable advantage for the behaviour study of the engine.

The multibody systems were considered for a long time as a group of rigid bodies linked with springs and dampers[10], [11], [23], and the structural flexibility of bodies are not considered. We should precise that the majority of these studied multibody systems are industrial robots, in general very massive and slow enough to justify the neglecting of the structural deformation. However, the changes of technology impose new constraints to multibody systems. The missions become more and more complex. In industry for example, we can remark an extensive use of fast and light arm robots, and in astronautics or aeronautics the flying objects possess long and light components. It is necessary in these cases to take into account the deformations of these components.

In aeronautics, the integration of the flexibility of the blades, particularly those of helicopters conduct to the lightening of the structures. This permits to these engines to be more operational for all possible missions. However, the growth of the rotating speed of the rotor increases the risks of vibration of these lights structures, and also the problem of sound nuisance. So the analysis of the dynamic behaviour of a helicopter which take into account the flexibility structural effect permits to construct numerical computational codes necessary to optimize the motion one part, and to establish a control strategy other part.

In this way, Houbolt and Brooks [9], conduct an original study, proposing an analytical model for the dynamic behaviour of an isolated blade in void. They establish the equations of motion of the blades without aerodynamic effects but with the conservation of the higher order terms due to cutting effort and to rotating inertia. This study opened the
way for several works based on analytic approaches. Vyas and Rao [22] improve this study by considering the variable speed of rotation of the blades. Culp and Murthy [5] subject an original resolving method of these equations by introducing an integration matrix.

Another technique of modelling of flying vehicles use the Kirchoff equations. This method were used firstly by Meirovitch [15], then Tuzcu and Meirovitch [21] to elaborate the dynamic model of robot with arborescent chains.

In this work, we establish the dynamical model of Quadrotor called XSF, which is assumed as a multibody system. We consider the XSF as a structure in cross form with a central body connected to flexible components.

For that, we consider firstly each substructure individually. We define kinematic quantity of each substructure in the local reference frame of central body. The mathematical equations traducing the motion of each one are obtained by using the Lagrangian variational method. The external loads acting on the quadrirotor are defined by the virtual work principle. We introduce the deformation of flexible components by a modal synthesis based on Rayleight-Ritz method. We retain uniquely the pure deflexion modes for the deformation of the flexible components. The aim object in this work is the development of a model combining lightness and accuracy with respect to computation time. The wished algorithm has to be precise, integrating the flexibility of components and the inertial coupling between the overall motion and the deformation. We conclude by numerical examples comparing our algorithm to a rigid body model of the XSF, to show how it is important to take into account the flexibility in such light flying objects.

## MATHEMATIC MODEL:

The kinematic description of all components of Quadrirotor will be established in the inertial reference. In order establish those mathematics expressions, one use the three orthonormal bases. One assume that the XSF Quadrirotor is a multibody system constituted by a central body ($B_c$) assimilated at a rigid cylinder which is solidly connected four flexible arms ($B_{f_j}$) with ($j=1,2,3,4$) assimilated each one at a deformable tubular beam. At extremity of each flexible arms is connected a rotor blades system ($B_{h_{jb}}$) assimilated at a embedded deformable thin rectangular beam free. A rotor blade system is constituted by two blades ($b=1,2$) symmetrically connected at a rotor. Consider the XSF shown in Fig.1. The kinematic of Quadrirotor is defined in

($\Re_g$) = {$T, \vec{x}_g, \vec{y}_g, \vec{z}_g$} assumed as inertial reference or Galilean reference. One defines by ($\Re_c$) = {$O_c, \vec{x}_c, \vec{y}_c, \vec{z}_c$} the local reference linked at central body, where $\vec{z}_c$ is the descending vertical axis, $\vec{x}_c$ the longitudinal axis belong of ($B_{f_j}$) for $j=2$ and $\vec{y}_c$ the longitudinal axis belong of ($B_{f_j}$) for $j=1$. In reality the presence of the flexible structural in the Quadrirotor system imply the superposition of flexible motion on the rigid motion for obtain the overall motion of each sub structure. A some formulation for dynamical analysis of the flexible structure use the floating reference approach for define the motion of the flexible body [4], [17]. In using the floating reference approach one can to define the boundaries conditions that permit to express the space functions of the modal deformation in using the discretization finite element method based on nodal deformed or the continuum modal function technique of the flexible body. Last ($\Re_{f_j}$) = {$O_{f_j}: \vec{x}_{f_j}, \vec{y}_{f_j}, \vec{z}_{f_j}$} denote the floating local reference attached at ($B_{f_j}$), $O_{f_j}$ is the contact point of each flexible arm with central body. ($\Re_{h_j}$) = {$O_{h_j}: \vec{x}_{h_j}, \vec{y}_{h_j}, \vec{z}_{h_j}$} denotes the floating reference of each propeller, $O_{h_j}$ the contact point with flexible arm. In this study we will not mention the case of the joints for each articulation treated in [13]. An other fact in the dynamic study of the mechanic system is take account of different attitude of central body. He is important for define the behaviour of Quadrirotor in flight. The rotation of central body around its centre inertia is define by Euler angle expressed by ($\varphi, \theta, \psi$) such as $\varphi$ is the rotation angle due to the roll motion, $\theta$ the rotation angle of the pitch motion and $\psi$ the rotation angle of the law motion. After to have define the sequence rotation such as: rolling-pitching and lace, on express by $A_c$ the matrix of passage of the local reference ($\Re_c$) to the reference inertial frame ($\Re_g$) given by:

$$A_c = \begin{pmatrix} c_{\theta}c_{\psi} & s_{\theta}c_{\psi} & -s_{\theta} \\ s_{\phi}s_{\theta}c_{\psi} - c_{\phi}s_{\psi} & s_{\phi}s_{\theta}s_{\psi} + c_{\phi}c_{\psi} & s_{\phi}c_{\theta} \\ c_{\phi}s_{\theta}c_{\psi} + s_{\phi}s_{\psi} & c_{\phi}s_{\theta}s_{\psi} - s_{\phi}c_{\psi} & c_{\phi}c_{\theta} \end{pmatrix} \quad (1)$$

$$c_{\varphi} = \cos \varphi; s_{\varphi} = \sin \varphi$$

$$\text{where: } c_{\theta} = \cos \theta; s_{\theta} = \sin \theta$$

$$c_{\psi} = \cos \psi; s_{\psi} = \sin \psi$$
# 1. Presentation of the XSF Model

<figure><img src="image_3.png" /><figcaption>Figure 1. Representation of the UAV.</figcaption></figure>

The XSF is a micro UAV (Unmanned Aerial Vehicle) represented by a quadrirotor of 68 cm x 68cm of total size. It is designed in a cross form and made of carbon fibre. Each tip of the cross has a rotor including an electric brushless motor, a speed controller and a two-blade ducted propeller. In the middle one can find a central cylinder enclosing electronics namely Inertial Measurement Unit, onboard processor, GPS, radio transmitter, cameras and ultrasound sensors, as well as the LI-POLY batteries.

The operating principle of the XSF can be presented thus:

Rotors {1} and {2} (see Fig. 1) turn clockwise, and the rotors {3} and {4} (see fig. 1) turn in the opposite direction to maintain the total equilibrium in yaw motion. The equilibrium of angular velocities of all rotors done, the UAV is either in stationary position, or moving vertically (changing altitude).

A characteristic of the XSF compared to the existing quadrirotors, is the swivelling of the support {6} of the rotors {2} and {4} around the pitching axis $\vec{x}_1$ thanks to two small servomotors. This permits a more stabilised horizontal flight and a suitable cornering (see Fig.2).

# 2. Dynamic model

The dynamic model of the XSF with rigid and flexible components is based on the Multibody Systems Dynamics (MBS). The modelling of rigid or flexible bodies in a MBS has been extensively used for robotics and terrestrial systems see [6] and [20]. In this approach we use the substructuration

methodology to study the dynamics of quadrirotor because of the presence of several elements. This substructuration method consists with the subdivision of the UAV in elementary bodies interconnected by kinematics joints such as shown in figure 2. In this study, we consider the central body as reference body for establishes the dynamic model of the whole system.

Let us consider the Fig. 2, the substructures {5-8} are a flexible arm, and the other substructures {5} and {7} are connected to the central body by a fix joint. However {6} and {8} is connected by a revolute joint around the $\vec{x}_1$ axis. The substructure {1-4} is the subsystem rotor-two blades, such as the rotor is the rigid body and the blades are a flexible body. At each free extremity of the arms {5}, {6}, {7} and {8} is connected the axis of the identical components {1-4} by revolute joints around the $\vec{z}_1$ axis. The whole of the system is in spangled form with a central body. We considered that the multibody system of the Quadrirotor helicopter is composed by a four identical systems made up a flexible fuselage and a rotor-two flexible blades system.

Thus, initially we begin our study with the study of system fuselage-rotor-blades. To establish the dynamic equations of this system. In this first step, we establish the mathematical equations of the fuselage and the rotor-blade system separately. In the second step we take account the revolute joint between the fuselage and the rotor-blades system by a multipliers Lagrange technical.

To establish the complete dynamic equations of the XSF we deduct the dynamical of the others bodies by symmetric properties such as that will be developed in the continuation.

## 2.1 Kinematics of MBS

<figure><img src="image_4.png" /><figcaption>Figure 2. Kinematic scheme of MBS.</figcaption></figure>
Let us consider that $O_c$ is the origin of the local reference frame of central body ($B_c$). The location of an arbitrary point $P_c$ of ($B_c$) with respect ($\mathfrak{R}_c$) is given by:

$$Y_c = R_{O_c} + d_c \qquad (2)$$

where: $R_{O_c}$ represent the position of origin $O_c$ and $d_c$ the local position of $P_c$. In using ($B_c$) as the reference body, one writes the vector position appertaining to other components of Quadrirotor such as:

$$\begin{align}
Y_{f_j} &= R_{O_c} + d_{cf_j} + u_{f_j} \tag{3} \\
Y_{h_{jb}} &= R_{O_c} + d_{cf_j} + d_{fh_j} + u_{fh_j} + A_{h_j}(d_{h_{jb}} + u_{h_{jb}}) \tag{4}
\end{align}$$

where: $u_{f_j}$ is the vector of flexible displacement of arm, $d_{cf_j}$ express the local position of the contact point between ($B_c$) and ($B_{f_j}$). $d_{fh_j}$ is the local position of the contact point between ($B_{f_j}$) and ($B_{h_{jb}}$), $u_{fh_j}$ its displacement vector due to the flexibility effect $d_{h_{jb}}$ express the vector of local position of blade and $u_{h_{jb}}$ its vector flexible displacement. Last use consider the position of blade express by Equ. 4, $A_{h_j}$ represent the matrix of cosine direction of the azimuthally angle of blade around of rotation axis given by:

$$A_{h_j} = A_{u_{f_j}} \cdot A_{\beta_j} \cdot A_{\psi_{h_j}} \quad \text{for } j=1,3 \qquad (5)$$

$$A_{h_j} = A_{u_{f_j}} \cdot A_{\psi_{h_j}} \quad \text{for } j=2,4 \qquad (6)$$

If we consider the rotation of the propeller with respect to its corresponding arm, we denote $A_{h_j}$ the rotation matrix of ($\mathfrak{R}_c$) with respect to ($\mathfrak{R}_{h_j}$). In fact such as shown in Fig. 2, the rotors {1} and {3} have the possibility of swivelling around the longitudinal axis of arm to which they are connected (i.e. around $\bar{x}_c$ local axis). That is permit to define the elementary rotation matrix $A_{\beta_j}$ ($j=1,3$) representing the swivelling angle $\beta_j$ such as: $-20^\circ \le \beta_{h_j} \le 20^\circ$.

$$A_{\beta_j} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & c_{\beta_j} & -s_{\beta_j} \\ 0 & s_{\beta_j} & c_{\beta_j} \end{pmatrix} \qquad (7)$$

with:

$$s_{\beta_j} = \sin \beta_j; \quad c_{\beta_j} = \cos \beta_j$$

Let us consider the hypothesis of small deformations for each deformable substructure. We consider uniquely the pure deflection motion of a beam. This permits to write the matrix $A_{u_{f_j}}$ expressing the rotation of element associate at the section due to the deformations of the beam given by:

$$A_{u_{f_j}} = I - \tilde{\mathcal{E}}_{u_{f_j}} \qquad (8)$$

where $I$ express the identity matrix

$$\tilde{\mathcal{E}}_{u_{f_j}} = \begin{pmatrix} 0 & u_{y_{f_j},x_{f_j}} & u_{z_{f_j},x_{f_j}} \\ u_{z_{f_j},x_{f_j}} & 0 & 0 \\ -u_{y_{f_j},x_{f_j}} & 0 & 0 \end{pmatrix} \qquad (9)$$

$u_{y_{f_j},x_{f_j}} = \frac{\partial u_{y_{f_j}}}{\partial x_{f_j}}, \quad u_{z_{f_j},x_{f_j}} = \frac{\partial u_{z_{f_j}}}{\partial x_{f_j}}$ correspond to the deformations of the beam in regard of $\bar{y}_{f_j}$ and $\bar{z}_{f_j}$ axis of ($\mathfrak{R}_{f_j}$).

$A_{\psi_{h_j}}$ is the matrix of direction cosines due to the rotation of the blades such as:

$$A_{\psi_{h_j}} = \begin{pmatrix} c_{\psi_{h_j}} & -s_{\psi_{h_j}} & 0 \\ s_{\psi_{h_j}} & c_{\psi_{h_j}} & 0 \\ 0 & 0 & 1 \end{pmatrix} \qquad (10)$$

with:

$$s_{\psi_{h_j}} = \sin \psi_{h_j} ; \quad c_{\psi_{h_j}} = \cos \psi_{h_j}$$

We can to write the velocity vector of a material point of substructure as follows:

$$\begin{align}
V_c &= V_{O_c} + \tilde{d}_c^T \omega_c \nonumber \\
V_{f_j} &= V_{O_c} + (\tilde{d}_{cf_j} + \tilde{u}_{f_j})^T \omega_c + \dot{u}_{f_j} \nonumber \\
V_{h_{jb}} &= V_{O_c} + (\tilde{d}_{fh_j} + \tilde{u}_{f_j})^T \omega_c + \dot{u}_{fh_j} \nonumber \\
&\quad + \overline{A_{h_j}(d_{h_{jb}} + u_{h_{jb}})^T} \omega_c + A_{h_j} \dot{u}_{h_{jb}} \nonumber \\
&\quad + \overline{A_{h_j}(d_{h_{jb}} + u_{h_{jb}})^T} \omega_{h_j} \tag{11}
\end{align}$$
where:

$\omega_c$ represents the angular velocity of local reference frame ($\mathbf{\Re}_c$) of ($B_c$) with respect to inertial reference frame ($\mathbf{\Re}_g$) which can to write as:

$$ \omega_c = G_c \Theta_c \qquad (12) $$

with:

$G_c$ the transformation matrix of angular velocity $\omega_c$ in derivation of Euler angles with respect to time such as:

$$ G_c = \begin{bmatrix} 1 & 0 & s_{\theta} \\ 0 & c_{\varphi} & -c_{\theta}s_{\varphi} \\ 0 & -s_{\varphi} & c_{\theta}c_{\varphi} \end{bmatrix} \qquad (13) $$

## 2.2 Discretization of the beam

We use the discretization based on the Rayleigh-Ritz method for define the vector of flexible displacement of each arm and each blade. We use the linear theory of beam such as the flexible displacement is only on the vertically and lateral direction i.e. $u_{ij}(x_{ij}, t)$:

$$ u_{ij}(x_{ij}, t) = \sum_{n=1}^{\mu} \phi_{t_{ij}}^{n}(x_{ij}) q_{ij}^{n}(t) \qquad (14) $$

$$ u_{ij}(x_{ij}, t) = \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \sum_{n=1}^{\mu} \begin{bmatrix} 0 \\ \phi_{y_{ij}}^{n}(x_{ij}) q_{y_{ij}}^{n}(t) \\ \phi_{z_{ij}}^{n}(x_{ij}) q_{z_{ij}}^{n}(t) \end{bmatrix} \qquad (15) $$

$\phi_{y_{ij}}^n$ and $\phi_{z_{ij}}^n$: represents the spatial modal functions define according to lateral and vertically axis represented by an interpolation function that will be developed in the next section.

$q_{y_{ij}}^n$ and $q_{z_{ij}}^n$: are the flexible coordinates systems.

In this work, we consider the flexible arm like a tubular flexible beam fixed with a mass with the end, the blades of the rotor are considered as a free fixed thin beam. The numerical values of the sharp modes selected are consigned in the following table such as:

<table><thead><tr><th>Mode retained</th><th>Flexible arm</th><th>Rotor propeller</th></tr></thead><tbody><tr><td>$k_1$</td><td>1.014</td><td>1.875</td></tr><tr><td>$k_2$</td><td>4.007</td><td>4.096</td></tr><tr><td>$k_3$</td><td>7.050</td><td>7.855</td></tr></tbody></table>

In considering the first modes, then we write the flexible displacement vector in compact form as:

$$ u_{ij}(x_{ij}, t) = \Phi_{ij} q_{ij} \qquad (16) $$

Where:

$$ \Phi_{ij} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ \phi_{y_{ij}}^{1} & 0 & \phi_{y_{ij}}^{2} & 0 & \phi_{y_{ij}}^{3} & 0 \\ 0 & \phi_{z_{ij}}^{1} & 0 & \phi_{z_{ij}}^{2} & 0 & \phi_{z_{ij}}^{3} \end{bmatrix} \qquad (17) $$

$$ q_{ij} = \begin{bmatrix} q_{y_{ij}}^{1} & q_{y_{ij}}^{2} & q_{y_{ij}}^{3} & q_{z_{ij}}^{1} & q_{z_{ij}}^{2} & q_{z_{ij}}^{3} \end{bmatrix}^{T} \qquad (18) $$

## 2.3 Kinetic energy

The kinetic energy of flexible body can be written as:

$$ T_{ij} = \frac{1}{2} \int_{(B_{ij})} (V_{ij})^T (V_{ij}) dm_{ij} \qquad (19) $$

$dm_{ij}$ is the elementary mass of body ($B_{ij}$) and $V_{ij}$ its velocity vector expressed above.

Let us consider the whole of system, the global expression of kinetic energy is given by:

$$ T = T_c + \left[ \sum_{j=1}^{4} T_{f_j} + \left( \sum_{b=1}^{2} T_{h_{jb}} \right) \right] \qquad (20) $$

Where:

$T_c$ represent the kinetic of central body ($B_c$), $T_{f_j}$ is the kinetic energy of each arm and $T_{h_{jb}}$ the kinetic energy of each blade.

The motion of multibody system is represented by a system of coordinates constituted of six freedom degrees of rigid motion introduced by a central body ($B_c$) and $6 \times (n_{f_j} + n_{h_{jb}})$ freedom degrees of

From the formulation of motion of the overall system depending of the variables of motion, one write the expression of kinetic energy in form:

$$ T = \frac{1}{2} \dot{r}^T M \dot{r} \qquad (21) $$

with:

$\dot{r}$ is the derivation of $r$ with respect to time such as:
$$r = \left[ R_{O_c} \quad \Theta_c \quad q_{f_j} \quad q_{h_{jb}} \right]^T \quad (22)$$

expressing the characteristics variables of motion
M is the mass matrix of system, which can be
formulated by the following expressions:

$$M = \begin{bmatrix} M_{RR} & M_{R\Theta} & M_{Rq_{f_j}} & M_{Rq_{h_{jb}}} \\ M_{\Theta\Theta} & M_{\Theta q_{f_j}} & M_{\Theta q_{h_{jb}}} & \\ \text{symmetric} & M_{q_{f_j} q_{f_j}} & M_{q_{f_j} q_{h_{jb}}} & M_{q_{h_{jb}} q_{h_{jb}}} \end{bmatrix} \quad (23)$$

The elements of mass matrix can be obtained from
equations (19) and (20) into take account the
formulation of variables of motion such as defined
in (22). One writes so that the expression of each
sub mass matrix in following form:

$$M_{RR} = \int_{(B_c)} dm_c + \sum_{j=1}^{4} \left[ \int_{(B_{f_j})} dm_{f_j} + \sum_{b=1}^{2} \left( \int_{(B_{h_{jb}})} dm_{h_{jb}} \right) \right]$$

$$M_{R\Theta} = \int_{(B_c)} \tilde{d}_c^T dm_c + \sum_{j=1}^{4} \int_{(B_{f_j})} \tilde{a}_{f_j}^T dm_{f_j} \\ + \sum_{j=1}^{4} \sum_{b=1}^{2} \int_{(B_{h_{jb}})} (\tilde{a}_{f_{nj}} + \tilde{a}_{h_{jb}})^T dm_{h_{jb}}$$

with:

$$\tilde{a}_{f_j} = \tilde{d}_{cf_j} + \tilde{u}_{f_j}$$

$$\tilde{a}_{fh_j} = \tilde{d}_{cf_j} + \tilde{d}_{fh_j} + \tilde{u}_{fh_j}$$

$$\tilde{a}_{h_{jb}} = \overline{A_{h_j} (d_{h_{jb}} + \tilde{u}_{h_{jb}})}$$

$$M_{\Theta\Theta} = \int_{(B_c)} \tilde{d}_c \tilde{d}_c^T dm_c + \sum_{j=1}^{4} \int_{(B_{f_j})} \tilde{a}_{f_j} \tilde{a}_{f_j}^T dm_{f_j} \\ + \sum_{j=1}^{4} \sum_{b=1}^{2} \left( \int_{(B_{h_{jb}})} (\tilde{a}_{fh_j} + \tilde{a}_{h_{jb}})(\tilde{a}_{fh_j} + \tilde{a}_{h_{jb}})^T dm_{h_{jb}} \right)$$

$$M_{Rq_{f_j}} = \int_{(B_{f_j})} \Phi_{f_j} dm_{f_j} + \sum_{b=1}^{2} \left( \int_{(B_{h_{jb}})} \Phi_{fh_j} dm_{h_{jb}} \right)$$

$$M_{\Theta q_{f_j}} = \int_{(B_{f_j})} \tilde{a}_{f_j} \Phi_{f_j} dm_{f_j} \\ + \sum_{b=1}^{2} \left( \int_{(B_{h_{jb}})} (\tilde{a}_{fh_j} + \tilde{a}_{h_{jb}}) \Phi_{fh_j} dm_{h_{jb}} \right)$$

$$M_{Rq_{h_{jb}}} = \int_{(B_{h_{jb}})} A_{h_j} \Phi_{h_{jb}} dm_{h_{jb}}$$

$$M_{\Theta q_{h_{jb}}} = \int_{(B_{h_{jb}})} (\tilde{a}_{fh_j} + \tilde{a}_{h_{jb}}) A_{h_j} \Phi_{h_{jb}} dm_{h_{jb}}$$

$$M_{q_{f_j} q_{f_j}} = \int_{(B_{f_j})} \Phi_{f_j}^T \Phi_{f_j} dm_{f_j} \\ + \sum_{b=1}^{2} \left( \int_{(B_{h_{jb}})} \Phi_{fh_j}^T \Phi_{fh_j} dm_{h_{jb}} \right)$$

$$M_{q_{f_j} q_{h_{jb}}} = \int_{(B_{h_{jb}})} \Phi_{fh_j}^T A_{h_j} \Phi_{h_{jb}} dm_{h_{jb}}$$

$$M_{q_{h_{jb}} q_{h_{jb}}} = \int_{(B_{h_{jb}})} \Phi_{h_{jb}}^T \Phi_{h_{jb}} dm_{h_{jb}} \quad (24)$$

Let us consider the different hypothesis on small
displacement of the deformations modular and the
length of the beam such as:

$$\frac{u}{x} \ll 1,$$

$$\Delta u \ll 1,$$

The elastic energy of each flexible component can
be expressed as:

$$U_{q_{f_j} q_{f_j}} = \frac{1}{2} q_{f_j}^T K_{f_j} q_{f_j} \quad (25)$$

$$U_{q_{h_{jb}} q_{h_{jb}}} = \frac{1}{2} q_{h_{jb}}^T K_{h_{jb}} q_{h_{jb}} \quad (26)$$

where:

$K_{f_j}$ and $K_{h_{jb}}$ represents the stiffness matrix of
each flexible arm and each blade respectively such
as:

$$K_{f_j} = E_{f_j} I_{f_j} \int_0^{l_f} \left( \frac{\partial^2 \Phi_{f_j}}{\partial x_{f_j}^2} \right)^T \left( \frac{\partial^2 \Phi_{f_j}}{\partial x_{f_j}^2} \right) dx_{f_j} \quad (27)$$

$$K_{h_{jb}} = E_{h_{jb}} I_{h_{jb}} \int_{0}^{l_{h_{jb}}} \left( \frac{\partial^2 \Phi_{h_{jb}}}{\partial x_{h_{jb}}^2} \right)^T \left( \frac{\partial^2 \Phi_{h_{jb}}}{\partial x_{h_{jb}}^2} \right) dx_{h_{jb}} \quad (28)$$

In considering the variables of motion of overall
system represented by the column vector *r*
expressed by Equ. (22), we can write the elastic
stiffness matrix of flexible component as follows:
$$ K = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & K_{f_j} & 0 \\ 0 & 0 & 0 & K_{h_{jb}} \end{bmatrix} \qquad (29) $$

## 2.4 Equation of Motion of the Body

From the Lagrangian formalism expressed by:

$$ \ell = T - U \qquad (30) $$

where:

$\ell$ is the Lagrangian of system, $T$ and $U$ the global kinetic and elastic energy respectively write above. In using the variational method based on the principle of virtual works, we can write the formulation of equations of motion in function of the coordinates system. We write the following equation:

$$ \frac{d}{dt} \left( \frac{\partial T}{\partial \dot{r}} \right) - \left( \frac{\partial T}{\partial r} \right) + \left( \frac{\partial U}{\partial r} \right) = Q^{ex} \qquad (31) $$

$Q^{ex}$: is the generalized force vector acting on the body that we will develop in the next section.

We obtain the dynamic equation in this form:

$$ \begin{aligned} & \frac{d}{dt} \left( \frac{\partial \ell}{\partial V_{o_c}} \right) + \tilde{\omega}_c \frac{\partial \ell}{\partial V_{o_c}} - A_c \frac{\partial \ell}{\partial R_{o_c}} = F \\ & \frac{d}{dt} \left( \frac{\partial \ell}{\partial \omega_c} \right) + \tilde{V}_c \frac{\partial \ell}{\partial V_{o_c}} + \tilde{\omega}_c \frac{\partial \ell}{\partial \omega_c} - (G_c^T)^{-1} \frac{\partial \ell}{\partial \Theta_c} = M \\ & \frac{\partial}{\partial t} \left( \frac{\partial \ell_{i_j}}{\partial \dot{q}_{i_j}} \right) - \left( \frac{\partial \ell_{i_j}}{\partial q_{i_j}} \right) + K_{q_{i_j}} q_{i_j} = F_{u_{i_j}} \quad (i=f,h) \end{aligned} \qquad (32) $$

where $F$ is the resultant force vector acting on the Quadrirotor and it includes the weight, the resultant aerodynamic force, obtained by integrating the aerodynamic density force over the entire surface of each blade; $M$ is the resultant moment vector about $A_c$ and it includes the moment due to aerodynamic forces only; $F_{u_{i_j}}$ is the force due to elastic phenomena of each flexible component of the Quadrirotor. Equations (32) are hybrid in the sense that the first two are ordinary differential equations, describing the translation and rotation of the Quadrirotor whole, and the last one is a partial differential equation, describing the elastic displacement of a typical point on each deformable component.

## 2.5 Generalized forces

The total force acting on the Quadrirotor in the central body reference frame is:

$$ F = mg + \sum_{j=1}^{4} A_{h_j} (F_{L_j} + F_{D_j}) \qquad (33) $$

where:

$F_{L_j} = \sum_{b=1}^{2} F_{L_{jb}} \cdot z_{L_{jb}}$: express the aerodynamics forces vector in the blade local reference frame.

$F_{D_j} = \sum_{b=1}^{2} F_{L_{jb}} \cdot x_{D_{jb}}$: express the drag force vector in the blade local reference frame

$g$: is the gravity acceleration vector

$m$: is the mass of Quadrirotor

$$ M = \sum_{b=1}^{2} \int_{(B_{hjb})} (\tilde{a}_{fh_j} + \tilde{a}_{hjb})^T (dF_{L_{jb}} + dF_{D_{jb}}) \qquad (34) $$

$$ F_{u_{f_j}} = \int_{(B_{hjb})} \Phi_{fh_{jb}}^T (dF_{L_{jb}} + dF_{D_{jb}}) \qquad (35) $$

$$ F_{u_{hjb}} = \int_{(B_{hjb})} \Phi_{fh_{jb}}^T A_{h_j}^T (dF_{L_{jb}} + dF_{D_{jb}}) \qquad (36) $$

We obtain the corresponding expression of generalized forces vector acting on the Quadrirotor helicopter

$$ (Q^{ex})^T = \left[ F \quad M \quad F_{u_{i_j}} \right]^T \qquad (37) $$

## 2.6 Aerodynamic forces and torques

In this part, we will define the characteristics of the aerodynamic forces and torques issued from the blade theory.

The blade behaves as a rotating wing. Each element of the blade dr. is in contact with the airflow with a speed $V_R$ and according to an angle of attack $\alpha$. One call pans the axisymmetric hooding of the hub, interdependent of the propeller in rotation. In the plan of the propeller, the pan is defined by the radius $S_p$.

Assuming in Fig. 3 $dL = dF_{L_{jb}}$ and $dD = dF_{D_{jb}}$.

Each elementary section of the blade of width $dr$ creates a lift $dF_{L_{jb}}$ and a drag $dF_{D_{jb}}$, such as:
$$ \left\{ \begin{aligned} dF_{L_{jb}} &= \frac{1}{2} \rho_{air} C_L V_R^2 dS \\ dF_{L_{jb}} &= \frac{1}{2} \rho_{air} C_D V_R^2 dS \end{aligned} \right. \qquad (38) $$

<figure><img src="image_5.png" /><figcaption>Figure 3. Description of forces applied on the blade.</figcaption></figure>

where $\rho_{air}$ is the air density, $C_L$ and $C_D$ represent adimensional coefficients of lift and drag depending mainly on the angle of attack $\alpha$. ($\alpha = \gamma - \Phi$).

The XSF will hover or move at low speed, we can assume that the thrust $F$ for a $b$ blades propeller can be written as:

$$ \left\{ \begin{aligned} F_{L_j} &= \frac{b}{2} \rho_{air} \omega_{h_j}^2 \int_0^{s_p} r^2 c(r) C_L(r) dr \\ F_{D_j} &= \frac{b}{2} \rho_{air} \omega_{h_j}^2 \int_0^{s_p} r^2 c(r) C_D(r) dr \end{aligned} \right. \qquad (39) $$

$$ \text{Or simply: } \begin{cases} F_{L_j} = k_L \omega_{h_j}^2 \\ F_{D_j} = k_D \omega_{h_j}^2 \end{cases} \qquad (40) $$

$k_L$ and $k_D$ are the coefficient of bearing pressure.

The computation of the lift coefficient is often complex. Moreover marginal swirls at the tip of the blade can modify significantly the theoretic values of some aerodynamic parameters.

It is thus essential to elaborate an experimental process, which will enable us to determine precisely the coefficient $k_L$ as well as the limits of validity of the equation (40).

<figure><img src="image_6.png" /><figcaption>Figure 4. Testing bench of the thrust.</figcaption></figure>

Let us denote $\overline{\overline{I}}_R \omega$ the kinetic moment of the rotor. The gyroscopic effects of the rotating elements are introduced in the model as follows:

$$ M_{Gy} = \overline{\overline{I}}_R \omega \wedge \Omega \qquad (41) $$

These gyroscopic effects appear especially if an air disturbance create a change in roll or pitch angle, and when we swivel the arm {5} and {6}. The gyroscopic moments are introduced in the global model (32) as generalised forces.

### 3. Simulation and numerical Results

The XSF is intended to move in an urban environment. Three kinds of manoeuvres are privileged: yaw rotation, vertical ascension, and translational displacement along the $\vec{x}_1$ axis.

In the first numerical test, we present the simulation of a vertical motion. Two models are considered. In the beginning we make the assumption that the XSF is a rigid body, subjected to external and gravity forces, and we compare the results with those of the full Lagrangian flexible model. The external forces applied on the XSF are gravity force and aerodynamic forces, which are function of the rotation speed of the rotor blades, here:

$$ \omega_{h_j} = 607.37 \, \text{rd/s}, $$

$$ M = 2.2 \, \text{kg} \text{ is the mass of the UAV,} $$

$$ l = 0.23 \, \text{m} \text{ the length of each arm.} $$

$$ E = 1.15 \times 10^9 \, \text{Pa} \text{ Young modulus} $$

$$ n_i = 2 : \text{shape mode retained for each flexible body} $$
The numerical simulation is based on the “semi explicit stable” Newmark method and developed in MATLAB®.

<figure><img src="image_7.png" /><figcaption>Figure 5. Comparison of aerodynamic forces acting on flexible and rigid model</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Comparison of global motion.</figcaption></figure>

The figure 5 represents the aerodynamic forces which acting on each rotor for the rigid and flexible model of XSF.

The figure 6 shows the global displacement of the centre mass of central body namely $O_c$. The blue-dashed line shows the influence of the flexibility in comparison with the same arm supposed rigid in a hand and with the rigid Newton-Euler model developed in [8] in the other hand.

In the second test we simulate a full manoeuvre of the XSF in an urban environment. In the beginning it should rise to an altitude of 2.55 m corresponding roughly to its cruise altitude when exploring villages. Then we swivel the rotors (1) and (3) to permit the horizontal displacement along the x-axis. After that we impose a “square” trajectory composed by four quarter of turn followed by

horizontal displacements. The different motions were controlled by a sample P.I.D. law.

<figure><img src="image_9.png" /><figcaption>Figure 7: Trajectory of the XSF in urban mission.</figcaption></figure>

The figure 7 shows the evolution of the position of the XSF centre of gravity in this mission. In a first approach we do not optimise the transition between one kinds of motion to another. This will be benefit to minimise the used energy.

## 4. Conclusion and future works

The equations of motion for manoeuvring flying structure in urban environment are non-linear due to the large body motion and the flexibility of some components. The model presented in this paper takes into account this non-linearity and includes the effects of flexibility, and the aerodynamic and gyroscopic effects. Then, the complete dynamic model of Quadrirotor is governed by six freedom degrees of motion of rigid body and the ($n_{f_j} + n_{h_{jb}}$) freedom degrees of motion of flexible components. That is an advantage in the numerical resolution of system.

We have considered one deformation mode for each flexible component. In this work, the flexible bodies were considered as a beam in pure deflexion. The dynamic model contain the non linear terms generated by the coupling between the rigid motion and the flexible motion. The contribution of flexible terms in the mass matrix conduct at no constant contrary to most models met in the literature.

In the part, consecrated to numerical simulation, we have showed the effect of flexibility of component on the dynamic behaviour of Quadrirotor such as shown by the figures. We have also showed the effect of flexibility in application of a control force in the rigid and flexible model case. The results obtained in the cadre of modelling by the relative coordinates method permit to observer the
effect of flexibility in this dynamic modelling of
XSF Quadrirotor. The model shows the influence of
the flexibility in reference to total rigid body models
widely used in this field. This model will be
completed later by the introduction of the
aeroelasticity at the blades.

Acknowledgement

Authors want to thanks the graduate students
Valentin Vala, Seneina Billel, Duaux Julien and Da
Silva David for their help in the design and
fabrication of the XSF, otherwise this project could
not be achieved.

References:

[1] J. A. Ambrosio. Computed-aided analysis of rigid and flexible mechanical systems, nato asi series. *Kluwer Academic Publishers Dordrecht*, 268, 1994.

[2] E. Altuge, J. P. Ostrowsky, and R.Mahony. Control of a quadrotor helicopter using visual feedback. *In Proceedings of the IEEE International Conference on Robotic and Automation*, 2002.

[3] F. Boyer, N. Glandais, and W. Khalil. Multibody dynamics based on a non linear euler bernoulli kinematics. *International Journal of Numerical Methods in Engineering*, 54: 27-59, 2002.

[4] F. Boyer and W. Khalil. Kinematic model of a multi-beam structure undergoing large elastic displacements and rotations, part two; kinematic model of an open chain. *Journal of Mechanism and Machine Theory*, 34: 223-242, 1999.

[5] J. D. Culp and V. R. Murthy. Free vibration analysis of branched blades by the integrating matrix method. *Journal of Sound and Vibration*, 155:303-315, 1992.

[6] M. Geradin and A. Cardona. A beam finite element non-linear theory with finite rotations. *Int. J. for Numerical Methods in Engineering*, 126,1988.

[7] R. J. Gibert. *Vibration des Structures*. Edition Eyrolles, 1988.

[8] E. J. Haug. Computer Aided Kinetics and Dynamics of Mechanical Systems. Allyn and Bacon, New York, 1989.

[9] J. C. Houbolt and G. W. Brooks. Differential equations of motion for combined flapwise bending, chordwise bending and torsion of twisted nonuniform rotor blades. Technical report, (-NACA-Repport-1346), 1958.

[10] W. W. Hooker and G. Margoulis. The dynamical attitude equations for an n-body satellite. *Journal Astronaut. Sci*, 12, 1965.

[11] W. W. Hooker. A set of dynamical attitude equations for an arbitrary n-body satellite having rotational degrees of freedom. *AIAA Journal*, 8, 1970.

[12] B. Jonker. A finite element dynamic analysis of spatial mechanisms with flexible links. Computer Methods in Applied Mechanics and Engineering, 76: 17-40, 1989.

[13] A. Mouhingou and N. Azouz. Dynamical analysis of a flexible multibody systems: Lagrangian approach. *In Proceedings ASME 8th Biennial Conference on Engineering Systems Design and Analysis*, Torino, July 2006.

[14] A. Mouhingou and N. Azouz. Modelling of a flexible quadrirotor. *In Proceedings ASME 30th Annual Mechanisms and Robotics Conferences*, Philadelphia, September 2006.

[15] L. Meirovitch. A unified theory for the flight dynamics and aeroelasticity of whole aircraft. *In Proceedings of the 11th VPI and SU Symposium Held Blacksburg*. SAE paper, 1997.

[16] M. Pascal. Vibrations analysis of flexible multibody systems. *Dynamics and Stability of Systems*, 6(3): 217-234, 1991.

[17] M. Pascal. Some open problems in dynamic analysis of flexible multibody systems. *Multibody System Dynamic*, 5: 315-334, 2001.

[18] W. Schiehlen. *Multibody Systems Handbook*. Springer-Verlag, Berlin, 1990.

[19] R. Schwertassek and S. Dombrowski. Modal representation of stress in flexible multibody simulation. *Non linear Dynamics*, 20: 381-399, 1999.

[20] A. A. Shabana. Computational Dynamics. Second Ed. Willey, New York, 2001.

[21] I. Tuzcu, *Dynamic and Control of Flexible Aircraft*. PhD thesis Virginia Polytechnic Institute and State University, 2001.

[22] N. S. Vyas and J. S. Rao. Equations of motion of a blade rotating with variable angular velocity. *Journal of Sound and Vibration*, 156: 327-336, 1992.

[23] J. Wittenburg. Dynamics of systems of rigid bodies. *Teubner*, Stuttgart, 1977.