Skip to content

Drag constant

In this section, you will experimentally determine the propeller drag constant \(k_d\).

Theoretical background

The propellers of a quadcopter act as aerodynamic surfaces, accelerating the airflow through them. This process consumes energy from the batteries and generates both thrust forces and drag torques on the quadcopter. As previously derived, the drag torque \({\color{var(--c2)}\tau}\) produced by a propeller is proportional to the square of its angular velocity \({\color{var(--c1)}\omega}\).

Drag Torque

\[ {\color{var(--c2)}\tau} = k_d {\color{var(--c1)}\omega}^2 \]

Where \(k_d\) is the drag constant (\(\text{N·m·s}^2\)).

Experimental procedure

You will measure the drag torque \({\color{var(--c2)}\tau}\) generated by the propellers using the yaw rotation rig, a test fixture that constrains every degree of freedom of the quadcopter(1) except rotation \({\color{var(--c1)}\psi}\) around its vertical axis.

  1. To secure the quadcopter, slide it in from the side and fasten it using two screws.

Device 2

By measuring the yaw angle \({\color{var(--c1)}\psi}\) as a function of time and knowing the quadcopter vertical moment of inertia \(I_{zz}\), which was derived in Exercise 3.2, it is possible to determine the total drag torque generated by the propellers through a dynamic analysis. To record the yaw angle over time, you should record the experiment using the slow-motion mode of your smartphone, with a stopwatch visible in the field of view.

Device 2 (Readings)

You must upload to the drone a program that spins propellers \(1\) and \(3\) at \(1000~\text{rad/s}\), and propellers \(2\) and \(4\) at \(2000~\text{rad/s}\). Since propellers \(1\) and \(3\) rotate clockwise and propellers \(2\) and \(4\) rotate counter-clockwise, a net torque will be produced, causing the drone to rotate in the counter-clockwise direction. For each quarter turn (\(90^{\circ}\)), record the elapsed time. Repeat the experiment three times and compute the average measurement.

To simplify the procedure, you may start and stop the propellers using the Take off and Land buttons in the Command Based Flight Control, located at the bottom-right corner of the Crazyflie Client.

Command Based Flight Control

Create a file named drag_constant.c inside the src/identification folder with the following code (1):

  1. Do not forget to update the motor model coefficients \(a_2\) and \(a_1\) (lines 8–9) with your estimation.
drag_constant.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#include "FreeRTOS.h"      // FreeRTOS core definitions (needed for task handling and timing)
#include "task.h"          // FreeRTOS task functions (e.g., vTaskDelay)
#include "supervisor.h"    // Functions to check flight status (e.g., supervisorIsArmed)
#include "commander.h"     // Access to commanded setpoints (e.g., commanderGetSetpoint)
#include "motors.h"        // Low-level motor control interface (e.g., motorsSetRatio)

// Motor coefficients of the quadratic model: PWM = a_2 * omega^2 + a_1 * omega
const float a_2 = 0.0f;
const float a_1 = 0.0f;

// Global variables to store the desired setpoint, the current state (not used here),
// and the computed PWM values for different motor speeds
setpoint_t setpoint;
state_t state;
float pwm_1, pwm_2;
float omega_1, omega_2; 

// Main application
void appMain(void *param)
{
    // Infinite loop (runs forever)
    while (true)
    {
        // Check if the drone is armed (i.e., ready to fly)
        if (supervisorIsArmed())
        {
            // Fetch the latest setpoint from the commander
            commanderGetSetpoint(&setpoint, &state);

            if ((setpoint.position.z) > 0)
            {
                // Set two different angular velocities for the motors
                // Motors M1 and M3 will spin at 2000 rad/s
                // Motors M2 and M4 will spin at 1000 rad/s
                // This configuration induces pure yaw rotation (spinning in place)
                omega_1 = 2000.0f;
                omega_2 = 1000.0f;

                // Convert angular velocities to PWM using the motor model
                pwm_1 = a_2 * omega_1 * omega_1 + a_1 * omega_1;
                pwm_2 = a_2 * omega_2 * omega_2 + a_1 * omega_2;
            }
            else
            {
                // If Z setpoint is not positive, apply minimal power to all motors (for idle spin)
                pwm_1 = 0.1f;
                pwm_2 = 0.1f;
            }
        }
        else
        {
            // If not armed, stop all motors
            pwm_1 = 0.0f;
            pwm_2 = 0.0f;
        }
        // Apply PWM to motors:
        // M1 and M3 get pwm_1 (corresponding to 2000 rad/s)
        // M2 and M4 get pwm_2 (corresponding to 1000 rad/s)
        // This asymmetric configuration results in yaw motion
        motorsSetRatio(MOTOR_M1, pwm_1 * UINT16_MAX);
        motorsSetRatio(MOTOR_M2, pwm_2 * UINT16_MAX);
        motorsSetRatio(MOTOR_M3, pwm_1 * UINT16_MAX);
        motorsSetRatio(MOTOR_M4, pwm_2 * UINT16_MAX);
        // Wait for 100 milliseconds before the next iteration (10 Hz control loop)
        vTaskDelay(pdMS_TO_TICKS(100));
    }
}

Follow the steps below to collect the experimental data:

  1. Ensure that the drone battery is fully charged
  2. Attach the drone to the fixture and place a stopwatch next to it
  3. Arm the drone by pressing the Arm button in the Crazyflie Client
  4. Start recording with your smartphone in slow-motion mode
  5. Turn on the motors using the Command-Based Flight Control
  6. Let the drone complete two full rotations and stop recording

After the experiment, extract the data and fill in the table below.

\({\color{var(--c1)}\psi}~(^{\circ})\) \({\color{var(--c3)}t}_1~(s)\) \({\color{var(--c3)}t}_2~(s)\) \({\color{var(--c3)}t}_3~(s)\)
\(0\)
\(90\)
\(180\)
\(270\)
\(360\)
\(450\)
\(540\)
\(630\)
\(720\)

Data analysis

Using the collected data, you should fit a curve relating the yaw angle \({\color{var(--c1)}\psi}\)(1) with the elapsed time \({\color{var(--c3)}t}\).

  1. Note that you must convert from degrees (\(^\circ\)) to radians (\(\text{rad}\)).

Drag Torque

Under the action of a constant torque, the angular displacement as a function of time is given by:

\[ {\color{var(--c1)}\psi} = \frac{{\color{var(--c2)}\tau_z}}{2 I_{zz}} \, {\color{var(--c3)}t}^2 \]

As already derived, this torque corresponds to the sum of the drag torques produced by each propeller \({\color{var(--c2)}\tau_i}\), which are proportional to the square of their angular velocities \({\color{var(--c1)}\omega_i}\):

\[ \begin{align*} {\color{var(--c2)}\tau_z} &= -{\color{var(--c2)}\tau_1} + {\color{var(--c2)}\tau_2} - {\color{var(--c2)}\tau_3} + {\color{var(--c2)}\tau_4} \\ {\color{var(--c2)}\tau_z} &= -k_d {\color{var(--c1)}\omega_1}^2 + k_d {\color{var(--c1)}\omega_1}^2 - k_d {\color{var(--c1)}\omega_3}^2 + k_d {\color{var(--c1)}\omega_4}^2 \end{align*} \]

Given that \({\color{var(--c1)}\omega_1} = {\color{var(--c1)}\omega_3} = 1000~\text{rad/s}\), \({\color{var(--c1)}\omega_1} = {\color{var(--c1)}\omega_4} = 2000~\text{rad/s}\), and \(I_{zz} = 4 \times 10^{-5}\,\text{kg.m}^2\), substituting into the previous equation yields:

\[ \begin{align*} {\color{var(--c1)}\psi} &= \frac{{\color{var(--c2)}\tau_z}}{2I_{zz}} {\color{var(--c3)}t}^2 \\ {\color{var(--c1)}\psi} &= \frac{\left(-k_d {\color{var(--c1)}\omega_1}^2 + k_d {\color{var(--c1)}\omega_1}^2 - k_d {\color{var(--c1)}\omega_3}^2 + k_d {\color{var(--c1)}\omega_4}^2\right)}{2I_{zz}} {\color{var(--c3)}t}^2 \\ {\color{var(--c1)}\psi} &= \frac{\left(-{\color{var(--c1)}\omega_1}^2 + {\color{var(--c1)}\omega_1}^2 - {\color{var(--c1)}\omega_3}^2 + {\color{var(--c1)}\omega_4}^2\right)k_d}{2I_{zz}} {\color{var(--c3)}t}^2 \\ {\color{var(--c1)}\psi} &= \frac{- \cancel{2}{\color{var(--c1)}\omega_1}^2 + \cancel{2} {\color{var(--c1)}\omega_1}^2}{\cancel{2}I_{zz}} k_d {\color{var(--c3)}t}^2 \\ {\color{var(--c1)}\psi} &= \frac{-1000^2 + 2000^2}{4 \times 10^{-5}} k_d {\color{var(--c3)}t}^2 \\ {\color{var(--c1)}\psi} &= \left(75 \times 10^{9}\right) k_d {\color{var(--c3)}t}^2 \end{align*} \]

Therefore, the most appropriate fitting model is a second-order polynomial with zero constant and linear terms. You should experimentally identify the drag constant \(k_d\) from this curve fit(1). Record the obtained value, as it will be used in subsequent activities.

  1. Tip: use MATLAB’s Curve Fitting Toolbox.

Results validation

Compare your experimentally obtained value with the one previously estimated with a ruler in Exercise 1.4. The drag constant \(k_d\) is expected to be on the order of magnitude of \(10^{-10}\,\text{N·s}^2\).