2.0 BACKGROUND

Reference 2 presents a review of available inertia measurement methods for various purposes. The main difference between the various methods is the use of either forced or free oscillations. Forced oscillation methods use an apparatus to force the test specimen into either translational or torsional oscillations^{(Reference de Jong and Mulder3)} with the inertial properties being measured indirectly by the force or moment required to move the test article. These methods can be very accurate^{(Reference Brancati, Russo and Savino4)} but require a complicated and expensive apparatus that is prohibitive for typical low-budget university research.

Free oscillation methods use some sort of translational or torsional pendulum. These will oscillate freely under the influence of gravity alone, once displaced from rest and released. Pendulum designs have been used with a single suspension wire as the simplest method or using multiple suspension wires to create either a translational or torsional pendulum^{(Reference Korr and Hyer5,Reference Previati, Gobbi and Mastinu6)} . The latter is more suitable for aircraft applications, because a multi-wire pendulum will hold the test article in a defined attitude.

The most widely used pendulum swing experiment for aircraft is the two-wire translational pendulum, used since the 1930s^{(Reference Miller7)}. An extension of the basic method^{(Reference Soule and Miller8)} discusses the aerodynamic effects affecting the pendulum motion due to the geometric properties of fixed-wing aircraft (buoyancy, entrapped air mass and the air mass that travels with the pendulum) and gives empirical methods to estimate the magnitude of those corrections, based on the geometric properties of typical aircraft components. These aerodynamic effects can add up to 20% on top of the true airframe inertias and are therefore very significant^{(Reference Soule and Miller8)}.

A further extension of the pendulum method^{(Reference Gracey9)} from 1948 indicates that a major source of difficulty of all previously mentioned methods is the requirement to accurately determine the pendulum length, which is the distance from the pivot to the pendulum CG. The difficulty lies with how the vertical CG position of the aircraft with respect to the pendulum pivot is to be determined. This method proposes to swing the pendulum with two different arbitrary lengths and use the two results to simultaneously solve for the pendulum lengths. This is potentially more accurate than measuring this quantity. A different method was developed in 1950^{(Reference Turner10)} to enable inertia measurements for large and heavy airframes. The pendulum method is impractical for such large airframes and it was replaced by a ground-based spring support. Otherwise, this method is similar to the pendulum method.

A novel method that uses multi-degree of freedom motion together with system identification^{(Reference de Jong and Mulder3)}, uses an apparatus which simultaneously allows rotation about the pitch and yaw axes and translation along the roll and pitch axes. The method does not require the measurement of the CG position of the test article, compared to the previously mentioned techniques. On the other hand, the required test apparatus is very complex.

There have been some existing references dealing specifically with small-scale UAVs^{(Reference Shakoori, Betin and Betin11-Reference Kotikalpudi, Taylor, Moreno, Pfifer and Balas15)}. Reference 11 compares different pendulum methods. The other references all use some form of the standard bifilar pendulum apparatus. On the very small-scale, Ref. 16 reports on a test of a small quad-rotor flight vehicle. The paper also performs a comparison with a CAD modelling method and concludes that considerable errors between physical experiments and CAD modelling are possible, unless extreme care in the modelling is taken. Another very promising, novel method, specifically for small-scale UAVs^{(Reference Bottasso, Leonello, Maffezzoli and Riccardi17)}, is based on a 3DoF pendulum, suspended from a three-axis gimbal. It uses system identification of the pendulum motion to estimate the entire inertia tensor at once. This method is the basis for the second technique to be presented in this paper.

As noted above, the results of the swinging tests of aircraft require substantial corrections due to aerodynamic effects. A very interesting report regarding these corrections looks into all kinds of error sources for a clock pendulum^{(Reference Lupton18)}. Some of these errors are very relevant for the current project and are indeed the same as reported by Soule and Miller^{(Reference Soule and Miller8)}. Other error sources, such as the change in gravitational acceleration due to to the position of the moon with respect to earth, are affecting only the long-term stability of a clock pendulum and do not need to be accounted for during the short experiment durations of the airframe swing tests.

Using system ID of the pendulum motion instead of timing the oscillation periods^{(Reference Soule and Miller8,Reference Gracey9)} has proven more reliable for the small inertias involved, because the system ID method uses the full motion data to estimate frequency and damping instead of characterising the motion from a few single data points from a timer. It is possible to obtain reasonably accurate data using a millisecond precision timer, but exploratory tests have shown that this requires many more experiment repeats than using the system ID method. Most recent implementations of inertia experiments use system ID for the data processing^{(Reference de Jong and Mulder3,Reference Mendes, van Kampen, Remes and Chu16,Reference Bottasso, Leonello, Maffezzoli and Riccardi17)} , and a similar approach has been taken for the presented work.

3.0 THEORY OF PENDULUM MOTION

3.1 1 DoF pendulum

The equations of motion for a rigid-body, 1DoF pendulum can be developed from the free-body diagram in Fig. 1. The test article with mass m, inertia I and frontal area S swings about pivot point O. The pendulum length or the distance between O and the CG of the test article is l. From Euler’s second law, the rotational equation of motion for the rigid-body pendulum can be written as

(1) $$\begin{equation} M=I_O\dot{\omega } \end{equation}$$

where M is the applied moment about O, I_O the inertia of the pendulum about O and $\dot{\omega }$ the rotational acceleration of the pendulum.

Figure 1. Free-body diagram of a physical pendulum.

For the system identification algorithm, the equations of motion need to be expressed in state space form, depending on the state derivative vector $\dot{\mathbf {X}}$ and the measurement vector Y. The state vector for the pendulum is X = [θ ω] ^T , where ω is the rotation rate of the pendulum and θ the attitude angle. Using Equation (1), the state rate equations for the rigid-body pendulum become

(2) $$\begin{equation} \mathbf {\dot{X}} = \left[ \begin{array}{c}\dot{\theta } \\ \dot{\omega } \end{array} \right] = \left[ \begin{array}{c}\omega \\ M/I_O \end{array} \right] \end{equation}$$

The applied moment is the sum of the moment due to the gravitational acceleration, which is driving the motion, and the opposing moments due to aerodynamic drag and bearing friction. The moment due to the gravitational acceleration about point O is

(3) $$\begin{equation} M_g = -mgl\sin \theta \end{equation}$$

The moment due to the drag about point O opposes the motion according to

(4) $$\begin{equation} M_D = - \bar{q}SC_D\times l \times sign(\omega ), \end{equation}$$\\

(5) $$\begin{equation} = -0.5 \rho S C_D \omega^{2} l^{3} sign(\omega ) \end{equation}$$

with the dynamic pressure $\bar{q}=0.5\rho V^2$ and V = ωl. The term sign(ω) ensures that the drag force is always opposing the direction of the motion. Friction in the bearings is small compared to the drag of the test article and will not be modelled separately. It has a similar effect on the motion as the drag; therefore, the identified drag coefficient will be slightly higher due to that friction.

So far, the equations of motion for the pendulum include the inertia of the pendulum about the pivot O. If the inertia of the test article about its CG is desired, the parallel axis theorem can be used to express the inertia I_O of the pendulum as

(6) $$\begin{equation} I_O=I+ml^2 \end{equation}$$

The measurement or output equations Y consist of the two states θ and ω that will be measured directly by the rig instrumentation. Using the above derivations, the final equations of motion can be written as

(7) $$\begin{equation} \mathbf{\dot{X}} = \left[\begin{array}{c} \dot{\theta } \\ \dot{\omega } \end{array} \right]=\left[\begin{array}{c} \omega \\ \frac{M_g+M_D}{I+ml^2} \end{array} \right], \end{equation}$$\\

(8) $$\begin{equation} \mathbf{Y} = \left[\begin{array}{c} \theta \\ \omega \end{array} \right] = \left[\begin{array}{c@{\quad}c} 1 0\\ 0 1 \end{array} \right] \left[\begin{array}{c} \theta \\ \omega \end{array} \right] \end{equation}$$

Inspecting Equation (7), some important implications for the experimental methodology can be immediately deducted. Firstly, because the expected inertias I of the test article are going to be small, it is paramount to make the pendulum length l as short as possible. For a small fixed-wing UAV, the parallel axis component of the total inertia of the pendulum, ml ², will always be large compared to the test article inertia I, and will dominate the magnitude of the denominator in Equation (7). If the pendulum is long, this will be even more significant. A long pendulum, therefore, makes it very difficult to identify the small inertia I with good precision. This problem has also been referred to in Refs 8 and 17.

Secondly, using knowledge of the system identification algorithm, it is also beneficial to estimate the full inertia of the pendulum about the pivot (I_O = I + ml ²) instead of estimating the inertia of the test article I in isolation. The parallel axis theorem can be applied after the estimate for I_O has be determined. The system ID algorithm perturbs each parameter by a small amount and calculates the significance of this perturbation onto the fit of the estimated system response to the measured data. Now, if Equation (7) is used in its stated form, a small perturbation of the parameter I will not have a significant effect on the magnitude of the denominator because ml ² is always much larger in magnitude. Perturbing I + ml ² instead gives a more robust response and the algorithm converges much faster and more accurately. Carefully considering those two issues during the experimental design is the first step to an accurate estimate of small-scale fixed-wing UAV inertial properties.

To determine the pendulum length l, the vertical CG location of the test article must be known precisely with respect to a reference point^{(Reference Miller7)}. This can be achieved by suspending the test article from a point longitudinally and vertically offset from the CG, as shown in Fig. 2(a). A vertical line from the suspension point to the longitudinal CG location will indicate the vertical CG position as shown in the figure. This can then be measured from any convenient reference line on the test article as indicated. The length l is then the distance from the pivot to this reference line plus the above measurement to the vertical CG of the test article. For a combination of objects, such as the test article mounted on a support frame, as shown in Fig. 2(b), the pendulum length l to the combined CG can be found from a moment balance about point O

(9) $$\begin{equation} l = \frac{m_Fl_F+m_{TA} l_{TA}}{m_F+m_{TA}}, \end{equation}$$

where the subscript F denotes the properties of the support frame and TA the properties of the test article.

Figure 2. Pendulum CG definitions.

The pendulum uses a support frame to hold the test article. This removes the need for any wire attachment points on the test article. Hence, after the swing test results are obtained, a final step is necessary to extract the inertial properties of the UAV about its CG. This step is the removal of the support frame inertias from the results. The support frame properties have to be determined beforehand by a separate experiment or other means, such as CAD modelling. A simple way of doing this would be to just subtract the inertias of the frame from the combined result, but this would lead to error because the inertia of the frame is typically measured about its own CG, which can be quite different from the combined CG location, as shown in Fig. 2(b). In addition, the CG location of the UAV is unlikely to be identical to the combined CG position as illustrated in the figure. To separate the inertial properties of the frame and the UAV accurately, it is therefore necessary to apply the parallel axis theorem, using the masses and dimensions given in Fig. 2(b). This leads to the equation for the test article inertia I_TA about its CG

(10) $$\begin{equation} I_{TA}=-[I_F+m_F\times l_2^2]-m_{TA}\times (-l_1)^2+I_{meas}, \end{equation}$$

where I_F is the frame inertia about its CG, m_F the frame mass, m_TA the test article mass, I_meas the combined inertia result from the system ID and the dimensions l_i as defined in Fig. 2(b).

3.2 3 DoF pendulum

The three-dimensional rigid-body pendulum, similar to that introduced in Ref. 17, is an extension of the 1DoF pendulum from the previous section. Figure 3 shows a sketch of the support frame with the axes definitions. The figure also contains an image of the 3DoF gimbal used to suspend the pendulum.

Figure 3. 3DoF Pendulum diagram and three-axis motion gimbal.

In Ref. 17, all equations were developed in the body frame of the pendulum with origin at the CG of the pendulum. As mentioned above, this is not ideal because it leads to numerical issues during the system identification. Hence all equations are developed in the earth fixed frame with origin at the gimbal pivot point O. This increases the stability of the system ID procedure. The body axes inertia tensor about the pendulum CG is computed using the parallel axes theorem as a second step.

The derivation method of the equations of motion for the 3DoF rigid-body pendulum is similar to the single-axis pendulum described before. Euler’s second law, now in three dimensions, describes the oscillatory motion of the rigid body as

(11) $$\begin{equation} \mathbf {M}=\mathbf {J}\mathbf {\dot{\omega }}+\mathbf {\omega } \times \mathbf {J}\mathbf {\omega }, \end{equation}$$

with M being the applied moment vector about the pivot point, J the inertia tensor of the pendulum about the pivot, and ω the vector of angular velocities p, q and r of the pendulum. For a typical aircraft with symmetry about the xz plane, the inertia tensor simplifies to

(12) $$\begin{equation} \mathbf {J}=\left[\begin{array}{c@{\quad}c@{\quad}c}I_x & 0 & I_{xz} \\ 0 & I_y & 0 \\ I_{xz} & 0 & I_z \end{array} \right] \end{equation}$$

The state space form of Equation (11) with states ω = [p   q   r] ^T , as required for the system ID algorithm, can be written as

(13) $$\begin{equation} \dot{\omega } = \mathbf {J^{-1}}(-\omega \times \mathbf {J}\omega )+ \mathbf {J^{-1}}\mathbf {M} \end{equation}$$

Similar to the 1DoF pendulum, the applied moment is the sum of the moment due to gravity and the moments caused by bearing friction and aerodynamic drag. The gravitational force vector acting at the CG is F _g = m g. The moment due to gravity can then be obtained by defining a vector R_{O, CG} from the pivot point to the pendulum CG and taking the cross-product

(14) $$\begin{equation} \mathbf {M}_G = \mathbf {R_{O,CG}}\times \mathbf {F}_G \end{equation}$$

The components of the vector R_{O, CG} in the earth fixed frame can be found by rotating the body axes vector R_{O, CG} = [0   0   l] ^T between the pivot point and the pendulum CG (as shown in Fig. 3) into the earth fixed axes using standard orthogonal transforms.

The damping terms are the moments due to bearing friction and the moment of the aerodynamic drag of the pendulum about the pivot point. The drag vector is the opposite of the velocity vector, which, due to the complicated motion of the 3DoF pendulum, constantly changes direction. Also, the drag will be different about each axis because the aircraft shape. Therefore, the drag is treated in component form for each axis, similarly to the bearing friction. No other aerodynamic force is expected to create a significant moment about the pivot point to influence the motion. Drag and bearing friction have a similar effect on the motion and can be combined into a single vector M_D = [M _{D, x}    M _{D, y}    M _{D, z} ], unless their separate numerical values are of interest. For this project, this was not the case, so the applied moment becomes

(15) $$\begin{equation} \mathbf {M} =\mathbf {R_{O,CG}}\times \mathbf {F} - \mathbf {M_D} \end{equation}$$

Expanding M into its components gives

(16) $$\begin{equation} \left[\begin{array}{c}M_x \\ M_y\\ M_z \end{array} \right] = l\left[\begin{array}{c}\cos \phi \sin \theta \\ -\sin \phi \\ \cos \phi \cos \theta \end{array} \right] \times mg \left[\begin{array}{c}0 \\ 0 \\ 1 \end{array} \right] - \left[\begin{array}{c}M_{D,x} \\ M_{D,y} \\ M_{D,z} \end{array} \right], \end{equation}$$\\

(17) $$\begin{equation} = mgl \left[\begin{array}{r}-\sin (\phi )-M_{D,X} \\ -\cos (\phi )\sin (\theta )-M_{D,Y} \\ -M_{D,Z} \end{array} \right] \end{equation}$$

The above equations require the attitude angles ϕ, θ and ψ of the pendulum with respect to the earth fixed frame. These states can be added to yield the final state vector

(18) $$\begin{equation} \mathbf {X} = [\phi \quad \theta \quad \psi \quad p \quad q\quad r]^T, \end{equation}$$

where the angles ϕ, θ and ψ are the attitude angles of the pendulum in the X, Y and Z axis, respectively.

The state rate equations for the Euler attitude angles can be found in any reference book on flight mechanics, such as Ref. 19. These, together with Equations (13) and (17) give the final equations of motion for the 3DoF rigid-body pendulum

(19) $$\begin{equation} \dot{\phi } = p + \tan (\theta )(q\sin (\phi )+r\cos (\phi )), \end{equation}$$\\

(20) $$\begin{equation} \dot{\theta } = q\cos (\phi ) - r\sin (\phi ), \end{equation}$$\\

(21) $$\begin{equation} \dot{\psi } = [q\sin (\phi )+r\cos (\phi )]/\cos (\theta ), \end{equation}$$\\

(22) $$\begin{equation} \dot{p} = I_{xz}\Gamma [q(I_x p+I_{xz}r)-I_y p q]+I_z \Gamma [q(I_{xz}p+I_z r)-I_y q r]-I_z M_x \Gamma +I_{xz} Mz\Gamma ,\qquad\quad \end{equation}$$\\

(23) $$\begin{equation} \dot{q} = [p(I_{xz} p+I_z r)-r(I_x p+I_{xz} r)]/Iy+M_y/Iy, \end{equation}$$\\

(24) $$\begin{equation} \dot{r} = -I_x\Gamma [q(I_xp+I_{xz}r)-I_y p q]-I_{xz}\Gamma [q(I_{xz}p+I_z r)-I_y q r]+I_{xz}M_x\Gamma -I_xM_z\Gamma ,\qquad\qquad \end{equation}$$

where Γ = 1/(I ² _xz − I_xI_z ). The measurement equation Y is

(25) $$\begin{equation} \mathbf {Y}= \mathbf {I}\times [\phi \quad \theta \quad \psi \quad p\quad q\quad r]^T, \end{equation}$$

where I is the identity matrix. The identified inertias about the pivot can then be transferred to the test article CG using the parallel axis theorem as before. The only terms affected are I_x and I_y , while for I_z and I_xz , the parallel axis theorem component is zero. The procedure of removing the support frame inertias from the solution is identical to the single-axis pendulum method.

4.0 AERODYNAMIC EFFECTS ON ESTIMATED INERTIAL PROPERTIES

During the swing experiments, the pendulum is immersed in a fluid (the surrounding air); therefore, the measured inertial properties will be affected by added mass due to the enclosed air, the buoyancy of the airframe in the surrounding air, and by the inertia of the air being accelerated by the pendulum^{(Reference Soule and Miller8,Reference Lupton18,Reference Brennen20)} . Hence, the measured inertia I_meas will be different from the test article inertia I_TA measured in a vacuum. Added mass is a general term for all aerodynamic effects that change the measured inertia of a body immersed in a fluid, as if the mass of the body itself was changed. The influence of buoyancy and enclosed air mass actually results in a change in weight of the test article, where buoyancy reduces the measured weight and the enclosed air adds to the measured weight of the test article, respectively. The added mass due to the inertia of the fluid accelerated by the pendulum motion is caused by a momentum change of the surrounding fluid due to an acceleration of the pendulum. It has the same form as a mass term^{(Reference Brennen20)}, hence the name.

As explained by Brennan^{(Reference Brennen20)}, a body moving through a fluid adds a certain amount of kinetic energy to the fluid. That kinetic energy T can be written for steady and rectilinear motion as

(26) $$\begin{equation} T=\frac{1}{2}m_{fluid}V^2 \end{equation}$$

or, if incompressible flow is assumed,

(27) $$\begin{equation} T=\frac{1}{2}\rho \zeta V^2 \quad \mathrm{where} \quad \zeta =\int _V \frac{v_x}{V}\frac{v_y}{V}\frac{v_z}{V}dV, \end{equation}$$

where V is the velocity of the body and the integral ζ is a measure of the volume of fluid affected by the motion of the body inside the entire fluid domain V. The resulting differences in velocity relative to V in the flow field are denoted v_i . The product ρζ is then the mass of fluid affected by the motion of the object. The integral ζ is constant for constant velocity V.

When the body accelerates or decelerates, as it constantly does during pendulum motion, the velocity V of the body changes and with that the kinetic energy T imparted on the fluid. This requires additional work to be done by the body, which is simply dT/dt. The rate of work done can then be expressed as F_DV, where F_D is an additional drag force. Assuming that ζ is constant, that is, the flow pattern does not change, the added drag, F_D , is

(28) $$\begin{equation} F_D=\frac{1}{V}\frac{dT}{dt}\quad =\quad \rho \zeta \frac{dV}{dt} \quad =\quad m\frac{dV}{dt}, \end{equation}$$

where the sign of the force depends on whether the body accelerates or decelerates. The added drag force F_D has the same form and sign as a force required to accelerate or decelerate the mass m of the body. Therefore, the term ρζ can be interpreted as an additional mass m_f of fluid that is being accelerated or decelerated by the body. This added mass m_f has an inertia about the axis of rotation and hence the inertia measured is

(29) $$\begin{equation} I_{meas}=I_{TA}+I_{m_f} \end{equation}$$

It should be observed that this added drag force is different from the ‘conventional’ drag force, which is proportional to the square of the velocity of the body. The added drag described here is proportional to the change in velocity of the body. Given the direct dependency on the fluid density ρ and the high density of water, this effect is very important especially for hydrodynamic problems^{(Reference Brennen20,Reference Lin and Liao21)} . However, it is also critical for the correct determination of airframe inertial properties as previously discovered^{(Reference Soule and Miller8)}.

For full-scale aircraft, the corrections due to the added-mass typically amount to 3% of the measured inertias due to buoyancy. The enclosed air adds around 5% in the X axis and is negligible in the Y axis. Another 20% are added in the X-axis due to I _mf , while errors in Y and Z, caused by I _mf , are about 5%^{(Reference Soule and Miller8)}. For small-scale UAVs, the enclosed air mass and buoyancy are negligible because the airframes only have small internal volumes and the volume of their structures is also small. Hence, the added mass due to these two effects will only add the equivalent of a few grams of weight to a 2–5kg airframe; therefore, these two corrections can be neglected. The third correction, I _mf , however, is very significant for small-scale UAV and thus requires careful consideration. For example, the correction due to I _mf in the X axis is 25.1% for the UAV used for this research.

Methods have previously been developed to estimate I _mf ^{(Reference Soule and Miller8,Reference Malvestuto and Gale22,Reference Gracey23)} . These methods are based on test data for flat plates and ellipsoids to model the airframe’s shape from these basic bodies. The corrections for the full test articles are then assumed to be the sum of the corrections for the separate components, ignoring potential interference effects between the parts. For full-scale aeroplanes, this appears to work quite well, with resulting inertia estimates within 2.5% or less of the true value^{(Reference Malvestuto and Gale22)}.

Applying the previously published methods to the given UAV geometry results in corrections that are an order of magnitude too small. It is unclear why this happens, since the datasets in Refs 22 and 23 are missing crucial information to reproduce the findings. Alternatively, Brennan^{(Reference Brennen20)} introduces a method for estimating the corrections based on results of a potential flow solver. Lin and Liao^{(Reference Lin and Liao21)} handled the same problem using modern fluid-structure interaction solvers. Both methods are conceptually very difficult and were beyond the scope of this project. Instead, geometrically similar models with known inertial properties, together with additional swing tests, will be used to determine the corrections due to I _mf for the given UAV.

It should be noted that none of the available publications on UAV inertial properties^{(Reference Shakoori, Betin and Betin11-Reference Kotikalpudi, Taylor, Moreno, Pfifer and Balas15)} acknowledge the requirement of these corrections. These reported results must, therefore, be treated with care.