QuadCopter Dynamics and Controls: Ramp-Up Reading…

I started to sketch a quadrotor model myself and derive the equations of motion, but I lost confidence and figured it would be best to sample the literature first. Below is a list of papers I found online. Click the links to see the papers.

All the papers are useful, but I found the, “PID vs LQ Control Techniques Applied to an Indoor Micro Quadrotor” by Bouabdallah et al. to be a most useful reference. It is becoming my primary basis for understanding and modelling the system I need to design a controller for.

PID vs LQ Control Techniques

modelling_of_quadrotor_minicraft

design_of_four_rotor_aerial_robot

Thesis KTH – Francesco Sabatino

State Space System Modeling of a Quad Copter UAV

Modelling of the ETH Helicopter

linearization

One of my goals is to develop a dynamic model and controller scheme as a platform for learning and applying, “Optimal Control and Estimation” techniques.

In reading these many papers I realized the multi-input model can be broken into a set of independent PID servos for, “level flight”. I decided I was more interested in understanding the coupled system in which roll, pitch, and yaw terms interact over a larger, “flight envelope” (larger roll and pitch angles) and not just small control action from level flight.

The Bouabdallah paper above doesn’t jump to model simplifications and linearizations. It maintains consideration for gyroscopic effects because the focus is on a lightweight, presumably highly dynamic platform. I’m also intrigued but the, “successive linearization” concept presented.

I’m going to get a better grasp on model linearization through this exersize. This means that I’m going to better understand how to modify the, “state transition matrix” (typically termed matrix ‘A’ for continuous control models and ‘F’ for discrete-time models). Hence the, “Linearization” paper above.

So I might be, “over-engineering” things here, for the sake of learning, but this is for fun so why not!

 




Dynamics and Controls Study Materials

I’ve enjoyed great opportunities to both study and work on mechanics, electronics, real-time control systems and supporting software for automated machines: Agricultural machine control for the most part, and some neat control systems problems in communications, sensors,  and other areas.

One thing I really like about, modelling a dynamic system is you need to be creative to arrive at an appropriate, “model”. It is just like a scale-model in that it isn’t the, “real thing”. It is something that represents the actual system, in a form simple enough to apply some math to, but not so complicated that you can’t solve the equations into some control electronics and software.

This, “modelling” step is my favorite part. It is the most important part. If it’s too complicated I can’t work with it. If it’s too simple, it means I’m probably ignoring something in the real world that would make any system design fail to do the job correctly. There really is not much math at this point, but a sketch and an accounting for what goes in, and what goes out: the stimuli you expect to make the system do it’s work, and also the forces of nature that want to move the system in ways you can’t always know, but need to manage somehow.

Because we’re dealing with, “dynamic” systems we’re talking about differential equations to represent them. This just means an equation that describes, given where something is (it’s current state), where it is going based upon some attributes and some forces acting on it right now.

It’s kind of like this: if you know where you are standing, you know how heavy you are, and how strong the kid is who is trying to push you out the door to drive him to a friend’s house, we might know how much closer to the door you will be one second from now.

These Differential equations are the mathematical model of a dynamic system. At this point of our model they might be too complicated for us to do much with, so we might need to decide if there are terms we can “ignore” because they will have a small effect, or maybe we need to “linearize” the system by making some assumptions about the range of action (motion) we expect relative to perhaps a typical operating point.

This part requires creativity too, and there are so many great examples of similar problems (especially now with the vast resources of the internet) that it is fun to, “borrow” techniques and sometimes combine techniques from different problems to a problem at-hand. Never forget, the best ideas are typically stolen. So, this, “modelling” of a system truly is like model-making when we think of plastic planes: we are creating a representation of a real thing, with many simplifications, only we are using math.

In December 2017 I decided to revisit system dynamics and control studies. I remember some of the old design “rules” and if I worked on some simple problems I could just grab some old notes, but I wanted a deeper (re)appreciation of some of the theory, and I wanted to strengthen my math skills again. In addition to the QuadCopter papers I shared in my last post, I’m studying the following texts.

Control Systems Basics: review of the old studies

This is introductory Controls textbook.

This is introductory Discrete-time controls text that builds on the text above. Most controllers today are implemented by a digital computer that samples sensors and outputs actuator commands on a, “sampling interval”.  We model this sampling and the math is a bit different (Z-transform in place of Laplace transform in the continuous time case). This is a good reference.

New Material…

I haven’t systematically studied “Optimal Control and Estimation” so after some old notes were challenging me I bought some books.

This looks like the classic text on this subject. I also found a great reference for working through the text:  Weatherwax_Stengel_Notes and some MATLAB samples  here

 

These last two were listed along with the Stengel book as some of the best, so I picked these up too. I’m a bit into this first one from 1970. I like some of the pre-personal-computer era books because they often tend towards math and hand calculations that promote deeper understanding.

 

 

 

 

I desired to brush-up on the Lagrangian method for deriving equations of motion. I found a link to a great reference online, and I bought a couple books from the guy! His name is David Morin from Harvard, and his classical mechanics materials look excellent…he also strikes me as a fun person in his presentation of materials.

 

I grabbed this one too…just for fun if I want to solve some problems. I’m not into Sudoku or crossword puzzles, so why not?

 

 

 

OK, enough with the books already. If you’re in school now or studied these topics before you may have simiar, suitable references.

 




QuadCopter DC Motor(BLDC) + Propeller + Gearbox Dynamics

I decided to start on the, “inner loops”. By this I mean the control loops for each propeller drive system. A Brushless DC Motor (BLDC) is standard for hobby-sized quadrotors, so let’s assume the motor type.

I haven’t sized anything yet: propellers or motors, but I can get the dynamic model figured out and put together a propeller speed control design.

Below I start with basic DC motor equations from just about any undergraduate dynamics textbook. I add the Bouabdallah model for the drag term, which for a simple motor would be the back-emf alone, but here I account for propeller and gearbox drag referred to the motor shaft as you can see below.

The motor-propeller speed model is non-linear. The Taylor Series expansion and evaluation about a nominal operating point is used to linearize it. As noted below, I’ll need a plan for sliding these constants and the compensator over the speed range. I can worry about that later.

I’m sharing this here because I want to illustrate more steps in the derivation of the final equation than the Bouabdallah paper has space to do. When I read papers often a final equation is given after it is stated that, “the following equation has been linearized”, for example. That is often a, “leap” for me. I don’t like to keep reading until I work it out for myself. Here I want to share each step.

I’ll leave this post to cover derivation of the propeller-motor equation with drive voltage as the control input.

BLDC




Timeout: QuadCopter Project…Why?

I didn’t do a very good job getting this, “QuadCopter” project kicked-off. I’ll probably interchange the term with, “quad-rotor” frequently. First, what is a “quad-rotor”? Well you probably know it’s a thing with four rotors. We buy it as a, “drone” and we fly it with a remote control or via a phone app acting like one. It is incredible what we can buy today for fun, and not too much money. Camera technology is awesome too.

I chose this gadget for a deep-dive because these things are popular. We can all remember seeing one fly, or maybe you own one. A quad-rotor is also a great example of a multi-output control system with multiple sensor inputs as feedback. We can go from simple to very complicated design ideas using this common vehicle as a platform for study.

How a Quad-Rotor Moves In the Sky

When you command a commercial quad-rotor drone to move up, down, right, or left relative to where it is in the sky, the machine processes that, “command” relative to the current location. Then it changes the speed of the propellers in such a way that the drone will move as you desire. A simple way of understanding it is as follows…

  • Two pairs of Propellers, opposite each other on the two arms that intersect at the center.
    • One pair rotates clockwise (CW – green above)
    • The other pair counter-clockwise (CCW – blue above)
  • Imagine driving one propeller of the green pair faster than the other green one.
    • if we left the blue pair equal, the body of the quad-rotor will rotate about a line between the two blue propellers above.
    • Same logic applies if we do the same to the blue pair, leaving the green equal: it will rotate about a line between the green propellers.
  • If we vary each pair more-or-less simultaneously we create a more complicated tipping action of the platform, but it’s basically just doing as described above to create a, “virtual” axis of rotation

One last motion: Yaw…

  • There’s one more less obvious motion: making the entire body we see above rotate about an axis coming out of the screen. This is, “Yaw” for this platform.
    • This relies on the counter-acting torque on the body of the quad-rotor as a result of the drive torque on the propellers: It’s Newton’s Third Law stating that for every action there is an equal and opposite reaction.
    • Here’s a way to visualize it for a single propeller. The plane wants to rotate the opposite direction. For a single-engine craft like this, the ailerons and rudder work to cancel this out.
    • The quad-rotor has four propellers. Here’s how we can make the quad-rotor, “spin” like a figure skater while remaining otherwise stationary
      • let’s say we start in a fixed spot in the sky, level, and not moving: hovering still.
      • If we speed-up the blue propellers above just a bit, and slow down the green just a bit so the total, “up” force is otherwise still the same…
        • the little “Extra” reaction torque on the body of the quad-rotor by the fast propellers compared to the slower propellers will cause it to spin one way.
        • If we flip it around and make the green a bit faster, the quad-rotor will, “spin” the other way.

No matter what the current roll, pitch, or yaw angle is at a given moment, if we increase or decrease every motor drive torque by the same relative amount the entire machine will move along the vertical axis of the drone, in theory.

Complications…

There’s a lot going on relative to non-linear propeller drive force equations, wind, and gyroscopic effects perhaps. The above explanation is simple enough for understanding the basics, but we rely on sensors to continuously measure and feeedback to the body controller so it can update the propeller speed controllers. This occurs very many times per second (probably around 50 times per-second would be a good rule for a small quad-rotor…3-to-5 times open-loop bandwidth, so conservatively say that is 10-20Hz times 3-5 is close to 50Hz…we’ll detail this when we get into picking a sampling rate later).

The presentation below does a good job of introducing the sensor concepts, and covers much of what I will likely repeat, but I am going to go into some more detail in a few areas to cover how to derive some of the equations, how to simplify the model to the degree shown here, or maybe not simplify it. For example, my recent post on the motor-propeller system modelling.

Also, this presentation looks to cover the simplified modelling and control using Proportional-Integral-Derivative (PID) controller pairs that act independent of one another to control pairs of propellers relative to, “level” flight.

That’s a fine and robust simplification for a stable camera-drone platform. I’m interested in learning what it would take to control one of these things over severe bank angles, or a larger, “flight envelope”. This is a good reference though, in addition to the papers I attached to my earlier post.


 




Quadrotor Propeller Speed Control: Commanding 4 props to stabilize the platform

In a recent post I shared a mathematical model for the combination brushless DC motor (BLDC), gearbox, and propeller model. A small quadrotor motor might not have a gearbox. It might have the propeller mounted directly to the shaft of a small hobby motor. A gearbox helps, “match” the motor to the, “load”: the propeller. Any gearbox is going to scale speed and torque, and add some inefficiency. We call these, “gain” terms: just multipliers in the chain of calculations. We can think about this again when we start shopping for actual motor-prop combinations for our (eventual) build. In the next post we will get (back) into our detailed equation for the motor-prop combination as we work on a propeller speed controller.

Back to basics…The speed of a propeller affects the exit velocity of the air behind it. The force the propeller offers in the forward direction is related to this fast moving air leaving the prop compared to the air in front of the prop. NASA has some great resources online. You can read about simple propeller dynamics here, from where I grabbed this image…

Jumping to conclusions, these equations tell us that the force offered by the propeller is proportional to the difference between the squared air-speeds across the propeller disk. The NASA link explains this is a simplified model, but we like it this way! Thank you NASA!

In the above equation for Force (in red) you can see that multiplier, “.5ρA”. That is the, “proportionality constant” that multiplies the squared air velocity difference to give us a force estimate. Propeller area is A (your basic area of a circle). The 0.5 comes from the dynamic pressure equation in white at the bottom of the black box above, and as the force equation is derived it flows through to the force equation.

The symbol ρ (rho) is the air density. It’s not constant. As an aircraft moves higher the density decreases, and as it moves through air of different temperatures the density changes. This is a simple model though, and our quad-rotor is not going to experience drastic air density changes such that the system could actually fail. This is a design consideration though: what are the limits of a system? For example, how high in altitude, into rarefied air (thin, not dense) could a propeller craft travel? At some point we simply could not produce enough force. Would a propeller get us very far in the thin air of Mars?

Back to earth, we are now going to get to the neat part of, “feedback” control: it let’s us get away with relatively simple physical models like the image above so long as we understand their limits and basic operation. We can be content with roughly knowing the proportionality constant above because we are going to use sensors on the quadrotor to, “feedback” tilts and position errors relative to where we want the platform to be.

Through some math these, “errors” will compute, “commands” for the 4 motors (armature voltage changes for each motor). These voltage changes will change the motor speeds. these speed changes will alter the squared-velocity terms in the above propeller force equations for the 4 propellers on the quadrotor. We don’t need to know the above equation precisely. In fact, the actual system will exhibit, “non-modeled dynamics” including, “non-linearities”: complicated physical phenomenon we don’t have clean-and-easy math for.

Instead, we basically ask for more-or-less motor speed many times per second on each of the four motors to keep it level, or to twist, or tip-and-move in one direction. We’re asking for relatively more-or-less, so we don’t need to be absolutely certain of all the parameters and values that govern the movement. There are things we can do wrong to make the commands oscillate and send the quad-rotor out-of-control. If we design our feedback control loops correctly we can guard against this.

This is all happening many times per second: sensors feeding-back to math that knows where we want to be, and more math to calculate individual motor commands based on the difference of where we are in the sky (and how we are tipped, tilted, or rotated) relative to where we would like to be, how level we want to be, and which direction we want our quadrotor to be facing. The, “math” will be simple arithmetic in a small computer on the quadrotor.

Remember the post with the dynamic model for the motor+gearbox+propeller? We ended-up with this differential equation for propeller speed:

Recall, “A” and, “C” are themselves functions of a nominal propeller speed, so the above equation is intended to apply around a certain propeller speed. The farther from that speed we drive the propeller by assuming pre-computed constant, “A” and, “C” the less reliable this equation will be. This means any controller we design on paper after assuming a nominal operating speed for the propellers could become invalid and cause our speed control to fail and the quadrotor to go out-of-control. We’ll need to determine how sensitive we are to this. If we need to worry about it we will modify our entire controller over the speed range.

We’ll think about this later. For now, the equation above appears below in block-diagram form. Notice the integrator on the right. When we implement this it will be as a, “difference equation” resulting from a Z-transform representation of what will be a sampled-data (computer controlled) system. For now I show the classic integrator there because you can see how the diagram implements the above equation.

We add a speed sensor block. You can see a reference input that will come from the main quadrotor flight controller: the amount of speed it wants out of the individual motor represented here. Again, this will be more-or-less speed error calculated many times per second. We need to design a, “compensator” that will produce armature voltage commands (“u”) based on the speed error; namely, the difference between the reference speed and what the prop speed sensor is telling us.

There will be four of these, “loops” running independently: one for each propeller on the quadrotor. The main flight controller will be supplying each with it’s own reference speed many times per second.

I’m going to end this post here. I started considering units and estimates for the parameters in the above equation (the A, B, and C terms) but getting into these details will be a good follow-up to this introductory post on the topic of controlling the motor speeds.

propellerSpeedControlLoop

 

 




Revisit the Motor+Gearbox+Propeller Model

 

It was fun to just mess with the motor-propeller differential equation and linearization in my earlier post. However, when I got into looking-up some rough parameter estimates and performing dimensional analysis (seeing if the units worked-out right for me in the equations) I ran into some problems.

I realized I was not only needing numbers, but I wasn’t comfortable with the units either! This means no controller design until I get a better understanding of the, “plant” (the thing we want to control: the propeller speed via the motor and gearbox if any).

Backup and learn more…

The picture for this post shows a Mr. Fred Ernest Weick (1899-1993). I show him because I learned of him as I dug into some details relating to this post. I got stuck trying to understand parameters in the motor+gearbox+propeller equation in my earlier post as I was moving ahead to a speed controller for each propeller.

Thanks to the internet I was able to find additional resources. Eventually I encountered a book by Mr. Wieck, and learned he was quite a force in the world of aviation. I found a PDF online of his 1926 book shown here. Mr. Weick did a lot of work on propellers and designed some aircraft. Safe civilian aircraft design appears to have been one of his motivations. I didn’t so much study his materials as get interested in his history. It is enjoyable to learn about aviation history.

Returning to our Motor+Gearbox+Propeller Equation

I derived this equation earlier to match what Bouabdallah presents in his paper (see his equation 14). I went step-by step to illustrate his result.

\frac{d\omega}{dt} = A\omega + Bu +C

Where…

A = -(\frac{K_m^2}{RJ} + \frac{2d\omega_0}{\eta r^3J})

B = \frac{K_m}{RJ}

C = \frac{d\omega_0^2}{\eta r^3J}

Trouble is, I had no idea where a couple of the terms originated. It was fine to work the equation earlier, but I did not fully understand the electro-mechanics through the gearbox to the propeller. I was, doing the math but not prepared to solve the problem!

I reworked the derivation many, many times after looking at more references. The biggest hang-up for me was the r^3  term in the equations. ‘r’ is the gearbox ratio. I could not figure out how it came to be cubed in these equations!

I kept re-reading, re-deriving, and looking for examples in other sources that would make it clear to me how this simple gearbox ratio gets cubed. In the end, it was very simple but it took me some hours off-and-on to finally, “get it”. I’ll share the process because this is one of my goals with these writings: share all the hang-ups and as many details as I can fill-in to illustrate not only the solutions to problems, but my struggles to understand all manner of details from the simple to the complicated. Sometimes they will be comical!

My, “sticking-point”

I was stumped by a small section of my primary reference paper for this project. In my earlier post I was happy to go through the linearization and get to the same, “answer” for the A, B, and C terms but I didn’t really understand one part of it: the propeller load on the motor through the gearbox.

Below is what I had to go on: steps 11-14 (on the left) in the paper together with the term descriptions to the right. The one step that reads, “By introducing the propeller and the gearbox models…” is what challenged my understanding. This is where \tau_d  in equation 12 becomes \frac{d}{\eta r^3J_t}\omega_0^2 in equation 13. I did not understand how this leap was made.

My understanding of the propeller load terms…

  • I was OK with the (d\cdot\omega_0^2) part that told me that the propeller put a load on the motor proportional to the square of the prop rotational speed.
    • This is consistent with the last post where the NASA equations had squared airspeed difference for total force out of the prop. It makes sense that we have a squared rotational term here for torque.
  • I had no idea where to find some example values for, ‘d’
    • I was not altogether confident I understood the context but I took it to imply the drag of the air resisting the rotation of a propeller.
  • How did the moment of inertia go from J to  J_t ?
    • I know when using motors with shafts, gears, and, “loads” we want to, “refer” all of the downstream complications to the motor shaft in one, “lumped” parameter so that our problem looks like a motor spinning a disk.
      • We, “lump” all the complications into an equivalent moment of inertia J_{equivalent}
        • I figured this is what J_t represents but was not sure how to arrive at it through the gearbox
  • Lumped inertia term aside, my main hang-up was the cubed gear-ratio in the denominator.
    • Typically we apply gear ratios as a fractional step-up in torque and down in rotational speed or vice-versa.
      • To see the gear ratio cubed as r^3 was a complete blind spot for me.

The Eye-opener

Reading another paper it suddenly clicked in my mind. It’s pretty simple. Most every time I am stuck I realize the resolution was simple later. This is why I don’t mind being stuck. I know if I just stay at a problem typically my mental block will clear.

Here the paper uses K_g instead of, ‘r’ so I’ll stick with the paper’s notation for a bit.

I knew these simple gearbox relationships but I had not systematically applied them to-and-from the motor to the propeller for the load model:

\omega_{gm}=K_g\cdot\omega_{gl}

\tau_{gl}=K_g\cdot\tau_{gm}

\tau_{gm}=\frac{1}{K_g}\tau_{gl}    tells us that the torque the motor sees is equal to the load torque at the propeller divided by the gear ratio (from the above gearbox equation).

This is the, “motor load” \tau_d term in equation (12) from our paper above, so here we see one gear ratio term in the denominator.

Now, on the propeller-side of the gearbox we have a simple relationship for torque as a function of propeller speed squared:

\tau_{gl}=d\cdot\omega_{gl}^2  BUT we also know that \omega_{gl} = \frac{\omega_{gm}}{K_g}

So it all simplifies to…

\tau_{d} = \tau_{gm} = \frac{1}{K_g}\cdot\tau_{gl}=\frac{1}{K_g}\cdot{d}\cdot\omega_{gl}^2 = \frac{d}{K_g}\cdot(\frac{\omega_{gm}}{K_g})^2 = \frac{d\cdot\omega_{gm}^2}{K_g^3}

 

\tau_{d} = \frac{d\cdot\omega_{gm}^2}{K_g^3}

There we have it! The K_g^3 term is our r^3 term! When we work this back through the equations we’ll get the A, B, and C terms the same as the Bouabdallah paper.

Here’s the paper that got me thinking correctly. This is an example of one method I use when I am stuck: just look at more-and-more source material until somebody’s treatment of a subject offers a different twist, or maybe includes an additional step, or diagram, or anything that helps you see it clearly. In the end, you can see I am combining the source materials. The Bouabdallah paper is still a primary reference, but this one and others are good support material. They help deepen understanding of each step in the process.

lab2_rotary_dynamics

 

Total moment of inertia as seen at the motor shaft

Let’s focus on line (3) in the PDF above so we can determine what our moment of inertia will be, referred to the motor shaft.

=\frac{J_l}{K_g}{\dot{\omega}}_{gl}    is simply replacing the \tau_{gl}  propeller-side torque with the J\frac{d^2\theta}{dt} definition of torque (the F=ma for rotations about an axis).

=\frac{J_l}{K_g^2}{\dot{\omega}}_{gm}   after we use the gearbox relationships above and replace the load-side \dot{\omega}}_{gl}  acceleration term with the motor-side equivalent through the gearbox.

That fractional term is the propeller load inertia as the motor, “sees” it on the motor side of the gearbox. We’ve referred it directly to the motor shaft. The total inertial load on the motor is then just the sum of the motor’s and this load inertia mapped back to the motor shaft.

J_{total}=J_m + \frac{J_l}{K_g^2} is then the inertia of the motor and the load referred to the motor shaft. In the case of the load, it is mapped-back to the motor shaft through the gearbox via the denominator term.

 

End Result of this exercise

My confusion over the r^3 term turned-out to be a blessing. I spent more time (re)working the electro-mechanical motor equations than I would have done. I accidentally became much more confident in all the parameters and my working knowledge of the derivation of the final differential equation and it’s terms:

\frac{d\omega}{dt} = A\omega + Bu +C

I really had TWO problems to overcome:

  1. My problem with the cubed gear ratio
  2. The total rotational inertia seen by the motor.

In my focus on settling #1 in my mind, I bumped into #2, and as indicated above I now have the equation I need to compute J_{total}

 

The, ‘d’ term in the propeller load equation

The trick with the, ‘d’ term is it needs to multiply an \omega^2 term to produce a torque. If I had studied any aeronautics I probably would have had a reference point for this, but I didn’t.

I eventually found a term in, “Modelling if the ETH Helicopter Laboratory Process” by Magnus Gafvert that appears to employ comparable motors and propellers for a small helicopter test stand. Here I found the following parameter listed:

“Aerodynamic torque for rotor R:    D = 2.91\times10^{-7} \frac{Nms^2}{rad^2}    (page 9, Gafvert) for what looks like a comparable hobby-sized propeller lab set-up. I can at least get started with this order-of-magnitude estimate. Those units will give me a torque (Nm) when I multiply by \omega_0^2 (nominal propeller speed around linearization point).

Gafvert references Mr. Weick in his paper, which is what sent me marveling at some engineering history. I could start my Aerodynamic studies with his 90 year-old book!

Conclusion

This post has been a great exercise in overcoming confusion. In my experience, grinding-through my own mental blocks and retreating from the blind alley is a big part of problem solving. Maybe that’s why we call them, “problems” and not immediately, “solutions”! I also come away with a much deeper understanding of the DC motor equations from my youthful days in school. I cemented good knowledge in preparing this report.




Propeller Blade, “Lift”

In  the last post I was able to clarify the load torque of the propeller applied through the gearbox to the motor shaft. I haven’t yet covered the main point of the propeller: LIFT! A propeller blade is just a spinning wing at each cross-sectional location.

The propeller force described in the NASA diagram in an earlier post summarizes the propeller-induced thrust force but it doesn’t tell us anything about the fluid mechanics around the prop blades.

Wings and Lift

My MIT Fluid Mechanics professor Ain Sonin would not be happy with me if I didn’t remember a very compelling demonstration he performed to illustrate, “streamline curvature” as the lift producer. The link on his name honors his passing in 2010. He was a young 72. I vividly recall his interesting lectures, impeccable sketches, approximation techniques, and his evident love of teaching. Thinking of what I learned from him, and what I should have better absorbed then, I want to cover, “lift” here.

Previously, we considered the torque load as the, “drag” on the motor. It makes sense that if we spin a propeller shaped like the example below the blunt nose and general form of the propeller requires some torque to spin it through the air. This is the drag you feel on your hand when you stick it out of your car window.

 

Let’s stick with the hand-out-the-window example. Whether we are actual kids or big kids, we all sometimes put our hand out the car window with our palm down. Then we tip our hand just a bit back and it feels like it wants to go up. Is this how a wing or a propeller (a spinning wing) creates lift? No!

A titled hand can rise with the force of air outside the car window but we cannot sustain flight (level-flight) or produce along-axis force from spinning a propeller unless the blades produce, “streamline curvature”! This is something hard to do with your hand out the car window unless you have a very unique shape to your hand and the angles are just right.

You are likely not producing, “lift” with your hand. you are just feeling the drag on an angle and that wants to push your hand up-and-back. If it weren’t attached to your arm it would just flip backwards. If an airplane wing attempted to climb by tipping itself backwards the plane would flip back like a mattress on top of a minivan going down the highway.

It’s about, “Streamline curvature” over the top surface

So what is really going on? In Professor Sonin’s class I was introduced to the fundamentals of Euler’s Equation, and its application in cylindrical coordinates along a, “streamline”. Streamlines are the blue lines flowing around the cross-section of a wing below. Bernoulli’s equation is what we learn to apply in basic fluid mechanics. It can be integrated from Euler’s equation. Euler gives us a deeper level of understanding.

If you look at the cross-sections of the propeller above you can see they look like a bunch of, “wing” cross-sections. The propeller is just a spinning wing. At any moment in time any of those propeller cross-sections above has streamlines flowing over-an-under similar to the simple wing diagram shown here. As it spins, part of the torque supplied by the motor goes into forcing these streamlines to curve around the blade. It is the nature of this curvature that gives us, “lift”!

For a wing, the, “weight” in the diagram is the weight of the entire plane, so the wings must produce lift to counter-act that weight, and extra lift to accelerate that weight upwards: to climb. When our propeller spins, we want it to pull a plane forward or a helicopter up (quad-copter in our case). A prop-driven airplane creates lift at the propellers to produce flight speed necessary to produce vertical lift by the wings! The prop, pulls the plane forward with, “lift” from the propellers and the plane stays aloft through, “lift” at the wings. the same principle is applied in two directions!

Our quadcopter will hover when the total lift from the four propellers matches the total weight. It will accelerate upwards or laterally when the total lift exceeds the total weight. The “force” we get out of a wing and a propeller is due to something you can see in the graphic above: the, “streamlines” on the top (B)  (the top surface of our quadrotor propellers) have more curvature than the streamlines on the bottom (C).

Euler’s Equation

The text below was in draft when I took Professor Sonin’s course at MIT. I still have a bound copy of the printout. I like this Text and I remember Sonin’s course for providing a more solid mathematical basis for fluid mechanics. Preparing for this post motivated me to (re)learn this material, so I bought a copy of Fay’s book for my shelf.

Recommended Reading

I’m going to gloss-over considerable detail (mostly because I would need to relearn too much right now!). I’m going to fill a couple gaps in the derivations that get us to, “streamline coordinates”. Otherwise I’m just organizing references here in hopes of creating a useful guide on this topic of lift as a result of, “streamline curvature”. I can’t improve the theory relative to these expert sources. I cannot come close to the genius of Euler and the father-son Bernoulli team! I hope I can make it more clear for we mere mortals is all.

Euler’s Equations

Euler’s Equations for inviscid flow will reveal what causes the lift once we get it into the right form. Inviscid flow is ideal, simplified flow assumed to slip past the wing with no surface interaction (no shear, so no boundary layer, eddies,  turbulence)…only the shape of the wing affects the flow, not the surface characteristics.

The Equations will reveal where the V^2 \Rightarrow\omega^2 for the propeller speed relationship to output force comes from. I realized I glossed over this topic in the  post with the NASA diagram where we saw the simplified equation with the V^2 term. In the last post where resolved some confusion around the motor+gearbox+propeller model for load on the motor I decided I needed to cover the entire point of spinning the propeller: “lift” from the blades, or aggregate thrust along the axis of the propeller.

The following PDF extracts key pages from Fay’s book (my old draft copy) that lead us to, “Euler’s Equation in Streamline Coordinates”. I retained a few pages of Bernoulli equation derivation. It is integrated from Euler’s equation along a streamline. It is a frequently applied formula. The typically overlooked perpendicular-to-streamline (normal) Euler equation tells us about lift from an airfoil!

James A. Fay: Euler's Equation

Example Application

This tutorial by the late Professor Sonin is excellent. I enjoy his neat sketches and remember them well, as I took notes in his class and attempted to match his detail and neatness.

Sonin: StreamLine Equations of Motion

It is not easy to digest the mathematical symbols, mutlivariable calculus, and terminology unless one has an excellent memory of their studies or practices these topics regularly. I will fill some gaps that troubled me as I re-derived some of the equations.

Cylindrical Coordinates

Although it might look daunting on it’s own, let’s accept Euler’s equation in cartesian (XYZ) coordinates. Then in Fay’s and Sonin’s material above we see we need cylindrical-streamline coordinates to get the equations of motion into some form that will enlighten us on this lift force we so desire.

Step-by-Step Derivation

The key to deriving cylindrical coordinates is to start with a representation of the unit vectors relative to a Cartesian system. I’m going to skip all the diagrams I should place here because you can look it up. I just want to offer a few observations on how elegantly simple the derivations are even thought they look complicated….

We start with radial and angular unit vectors in an X-Y plane, and we know Z in the Cartesian system is same as Z in Cylindrical

i_r = cos(\theta(t))i + sin(\theta(t))j

i_{\theta} = -sin(\theta(t))i + cos(\theta(t))j

i_z = i_z

A big difference between cylindrical coordinates and Cartesian coordinates is the rotation of the ‘r’ and ‘\theta‘ axes in time, as our, “fluid particle” moves in \theta. Hence showing \theta as a function of time in the above equation.

We’re going to need derivatives of the unit vectors, and here’s what they are:

\dot{i_r} = (-sin(\theta(t))i + cos(\theta(t))j)\dot\theta  but look at the trig in parenthesis. It’s our term for i_{\theta} above so…

\dot{i_r} = \dot\theta i_\theta

And

\dot{i_\theta} = (-cos(\theta(t))i -sin(\theta(t))j)\dot\theta   but again we see a familiar expression from above so this simplifies to…

\dot{i_\theta} = -\dot\theta i_r

It is these moving axes that give us the, “extra” terms in the cylindrical equations of motion. We get to them simply by taking derivatives of position and velocity.

We start with the position of a particle in cylindrical coordinates.

P = r\cdot i_r + z\cdot i_z         By chain rule we differentiate position to velocity as…

V = \frac{dP}{dt} = \dot r\cdot i_r + r\cdot \dot i_r + \dot z \cdot i_z + z\cdot \dot i_z

But we have terms for the derivatives of the unit vectors from above, and \dot i_z = 0 so velocity simplifies to…

V = \dot r\cdot i_r + r\dot \theta i_\theta + \dot z \cdot i_z

And because I don’t want to skip steps we differentiate Velocity to acceleration again by the simple chain rule of differentiation…

a = \frac{dV}{dt} = \ddot r\cdot i_r + \dot r\cdot \dot i_r + \dot r\dot \theta i_\theta + r\ddot \theta i_\theta + r\dot \theta \dot i_\theta + \ddot z \cdot i_z + z \cdot \dot i_z

and we again replace those unit vector derivative terms and \dot i_z = 0 to get, with a bit of rearrangement…

a = \frac{dV}{dt} = (\ddot r - r\dot \theta^2) \cdot i_r + (2\dot r \dot \theta+r\ddot \theta) \cdot i_\theta+\ddot z \cdot i_z|

Remarks on Cylindrical Coordinate Derivation

I must admit that unless I derive these equations as I did here I can fall prey to thinking there is something, “magic” about all those terms in the cylindrical coordinate system compared with our more intuitive Cartesian (XYZ) system. Perhaps for me more than for you I captured the derivation above so I can remind myself there is no magic to it, but simply two rotating unit vectors and the chain-rule of differentiation.

Streamline Coordinates

These are just a slight twist on Cylindrical coordinates. We relate cylindrical coordinates to streamline coordinates as follows…

i_n = -i_r       The perpendicular-to-streamline or, “normal” to flow unit vector.

i_s = i_\theta          The along-path unit vector, tangent to streamline unit vector.

i_b = i_z          By right-hand-rule, normal to the above two unit vectors: The, “up” axis of a cylindrical system. Out-of-page axis here. No fluid dynamics along this axis.

\frac{\delta}{\delta n} = -\frac{\delta}{\delta r}    As the above relationship between streamline normal and cylindrical radial unit-vectors, this is the derivative relationship.

\frac{\delta}{\delta s} = \frac{\delta}{r\delta \theta}      Along-streamline derivative relationship

That last denominator might be confusing but it’s just as simple as spanning a distance ds with angular change d\theta operating on a radius r.

If you now go back to the Fay PDF and accept Euler’s Equation in Cartesian coordinates on page 90-91 you can use the steps above to get confident with the cylindrical acceleration representation, also on page 90. From there you use the relationships just above to adjust the cylindrical representation to the streamline representation and you get to Fay’s page 105 below. We’re getting close!

Understanding, “lift” from Streamline Coordinates

The key equation for us is

-\frac{\delta p}{\delta n} = \rho\frac{V(n)^2}{R(n)}

I added the emphasis on V and R being likely functions of n, so integrating to estimate a pressure difference would be non-trivial. Sonin’s paper gives us a chance to solve the n-direction equation with models for Velocity and streamline radius as a function of n. See the paper to attempt this derivation for flow over a simple hill. I have not yet been successful arriving at equation (25) by using (23) and (24) to solve the differential equation (21) (I think it’s an integration-by-parts trick I am failing on).

As Fay and Sonin state, this is not a practical method for actually calculating a lift force. It is over-simplified, and we seldom have representations for V and R as a function of n over a range of s-vectors tangent to the flow over a wing. If we did, we’d still need to integrate them. However, this treatment tunes us into the physics of the problem. That’s our goal here: physical understanding, even if exact solutions elude us.

The approximation below from my 1995 Sonin lecture notes show us what we need to know in order to appreciate why our prop-speed-squared (\omega^2) is our lift producer. The notes were taken here, which I remember because I recall daydreaming as I watched a guy mow the grass outside the lecture hall window one fine autumn day.

Lift lecture notes (must not have been staring out the window at this time)

Understanding Approximations for Lift

The n-direction equation tells us what we need to know to understand lift generated by airfoils (thrust from propellers: spinning airfoils). As Sonin’s paper above states, “The n-direction equation states that when there is flow and the streamlines curve, the sum p+\rho gz (which is constant when the fluid is static) increases in the n-direction, that is, as one moves away from the local center of curvature.”

If we imagine integrating from the surface of the wing to far above the wing the streamlines way out at our approximated infinity would have a very large radius while the airspeed is V_\infty

If we just assume airspeed is V_\infty everywhere (a crude approximation because we are speeding-it-up over the wing) but let the streamline radius get larger from top-surface radius R as we move away from the wing then we can estimate a low pressure on the top of the wing proportional to \rho V_\infty^2 as indicated in the lecture notes above.

Remember that if our propeller is spinning with angular velocity \omega. At any point along the propeller where we could take a cross-section to sketch an airfoil diagram the “V” is just propeller rotational rate times the radius from the propeller axis to the location of our cross-section.

V_{cross-section} = \omega_{prop}\cdot r_{cross-section}

At this particular propeller cross-section we could estimate the free-stream velocity-squared as

V_\infty^2 = V_{cross-section}^2 = (\omega_{prop}\cdot r_{cross-section})^2 =\omega_{prop}^2\cdot r_{cross-section}^2

There’s our propeller \omega^2 term!

We now see how propeller rotational speed squared (\omega^2) is the variable responsible for, “lift” generated by a propeller blade. It is proportional to a simplified, assumed freestream velocity V^2 at any point along the blade where we take a cross-section as above. This is a highly simplified model, but it gets us to the essence of an airfoil.

Observe this on your next flight

Next time you land or take-off in an airplane observe the flaps that extend down-and-back from the rear of the wings. They create a shorter wing curvature radius and a larger wing surface. This increases lift. It is required because the plane is going much slower at take-off and landing than when cruising so the pressure difference across the same wing would be less at low speed, and it would be difficult to keep the plane, “up”.

So the modifiable airfoil extends the flaps, creates a larger surface over which the pressure difference acts, and puts more curvature into the airfoil to increase the pressure difference at low airspeed. This occurs at the expense of drag however. Given it is only needed at low speeds the flaps are retracted for cruising, where V_\infty is higher so we don’t need the curvature and extra area the flaps provide.

There remains that \rho term for air density, and we know it is higher in the troposphere down near the ground than up at 30,000 feet and higher. We’re going to ignore this though because our quad rotor is going to operate near the ground. Planes are more sensitive to, “stall” (failure of the lift from wings) at high altitude due to the rarefied (low density) air. Airspeed and angle-of-attack are critical factors, “up there”.

Conclusion

We now know that when our quad-rotor main controller asks our 4 motors for more-or-less speed many times per second it is the speed-squared (\omega^2) that gives us lift! I know I wrote this before, and have taken some deeper-dive detours, but from here we should get back to the prop speed control loops and soon onto the quadrotor main platform stabilization control.

This post has been a great journey back to interesting topics I once studied. Just as through the last post, the opportunity to produce this write-up has forced me to (re)learn to more depth and understanding than I learned it the first time! I hope the references above and any gap-filling explanations I offer help convey the, “lift” phenomenon!

 

 




Quadcoptor Platform Equations of Motion: Dynamic Model

It’s been over one month since I last posted. Back then I took a journey into fluid mechanics to offer a glimpse into how propellers work. I was motivated to do this as I was previously working-out the electric motor and propeller dynamics when I started to recall my mechanics studies.

It’s a good time to revisit my method for going at this, “system” problem of quadcoptor design.

Quadcoptor Project Review:

Basic Design Objectives

I added an objective to try a design that doesn’t assume the quadcoptor will always be level.

  • Have you noticed that camera drones don’t bank? There are some simplifications to the design if you want to maintain a horizontal platform at all times.
  • I’m curious to model a system that can operate away from the horizontal plane, at least in theory. I’m designing for a larger, “flight envelope”.

Sub-System Breakdown

  • When we get to a physical build there will be electrical, mechanical, and software sub-systems but at this point I’m looking at the system on paper:
  • How does the problem break-down into clear input-output relationships? Think of these as layers:
    • One level operates with certain data and makes a decision to command another, lower layer.
    • That lower layer operates with the above command and other information (sensors) to command what is it responsible for.
    • …and so on producing…
    • Physical action of the quadcoptor that is sensed again by each layer in this process which repeats itself.

This repeating is what makes these layers, “loops”. They process inputs, command outputs, sense the environment as inputs…and go again in a, “loop”. The sampling-rates are the time between successive trips through each loop at each layer.

Layers or, “Loops” I am working within

  1. “Flight Planning”: I’m not even into this yet so ignore it for now:
    1. Whether we use a joystick or upload some coordinates our craft will follow based on GPS.
  2. “Platform Control”
    1. The Quadcopter flight dynamics: the complete system in the sky that can move up, down, sideways, and rotate.
    2. Perhaps bank and perform some complicated actions predictably. This will be a, “flight envelope” or range-of-motion that likely exceeds a commercial, “camera drone” and for which we can come to understand the mechanics and the limitations.
    3. Sensor inputs to the platform. This will involve sensors for position and orientation: we’re not there yet but it will involve measuring and calculating positions and angles we’re getting into now.
    4. Outputs from the, “platform controller”: Propulsion
  3. “Propulsion”: I committed to a, “quad-copter” so we’ll have four propulsion units
    1.  Kinematics: how the placement and geometry of the 4 propulsion units affect the platform with their output thrust.
    2. Kinetics: producing the thrust from each output unit.

Earlier Posts Covered Propulsion System Modeling

If you look back you’ll see we covered motor-propeller electro-mechanics including a foray into fluid mechanics of lift from propellers. We now know the thrust producer is going to be a function of propeller speed. Our, “platform controller” is going to calculate how it wants to move, determine what combinations of the four propeller speed changes it needs to make to attempt this movement. Then it will command each motor-propeller propulsion unit accordingly. This, “loop” will likely run at near 10Hz or-so (this will be explained in a future post).

These four propulsion units each have a speed control loop to maintain the commanded propeller speed. These are the innermost, “loops” of control. This loop will run faster than the platform control loop (some tens-of-Hz. Again, we’ll cover this later). The specifics will depend on the motor speed controller we buy. We will not likely make our own motor speed controller.

No motors, propellers, sensors, or controllers for each propulsion unit have been selected. We understand the electro-mechanics though. We’re not going to shop for any parts until we get all the way through the system design. We’re headed there now!

Quadcoptor Platform Control

In my reading and on paper I’ve been studying in this area from the beginning even though I’ve only been writing about the propulsion unit. System design is an iterative process. In the mind it is all mixed-up and everything is interrelated.

Quadcoptor Equations of Motion

Degrees-of-Freedom

Our quad operates with 6 degrees-of-freedom (DOF). Its center-of-mass (historical term is center-of-gravity CG so this will be our abbreviation) moves up, down, and sideways relative to an “inertial reference frame”. It also tilts-and-twists as described by three angles below. That term, “Inertial reference frame” just means fixed relative to the environment our object will be moving: not accelerating.

A typical choice of reference frame for a hobby-grade quadcoptor is the, “take-off” point as the zero-point reference, an X-axis pointing north, Y-axis pointing east, and Z-axis pointing down into the earth via the right-hand-rule relative to X and Y axis. Point your right arm north with palm facing east. Curl your fingers towards the east. That is the right-hand-rule from X-to-Y axis. Extend your thumb. It points down. That’s your local, “Z” axis. Let’s call this our North-East-Down (NED) reference frame.

Reference Frame Assumptions

Our craft will operate within sight of this point. We will assume a flat earth. As our quadcoptor moves in the sky the CG position will be determined as an (X,Y,Z) coordinate relative to the inertial reference frame we chose to pin to the ground: this NED frame. Technically, this point on the earth is accelerating as the earth rotates but we’re assuming the entire local flight space is rotating through the day as we all do here on earth, and we accept this as a fair, “inertial frame”.

The quadcoptor is not a simple point mass. If you hold a ball-bearing and move it around you can think of it as a point-mass, and even though you know a bearing twists and rolls ignore that and assume it doesn’t rotate at all. The quadcoptor ball-bearing equivalent places all the mass of the object at the CG for the purposes of the above (X, Y, Z) position relative the the inertial NED reference frame. However, the tilts and rotations of the quadcoptor define the, “orientation” of the entire assembly in the sky.

For aircraft there are three angular rotations commonly used to describe orientation:

\phi = roll

\theta = pitch

\psi = yaw

These gives us our remaining 3 degrees-of-freedom, for a total of 6 relative to the earth-fixed NED reference frame:

(X,\: Y,\: Z,\: \phi,\: \theta,\: \psi)

The, “Body Frame”

Figure 1-2 near the end of the Blakelock Chapter 1 extract below offers a helpful illustration and legend describing the body frame, fixed NED frame, and Euler angles. Here’s the figure. See the PDF below for a detailed description on pages 15 and 16.

The sketch below illustrates a coordinate system fixed to the body of the quadcoptor. The X-Axis is defined to point always from the CG to the location of the first propeller. The Y-axis points towards what I label as the fourth propeller axis, and again by right-hand-rule Z is orthogonal to body X-Y plane by the right-hand-rule.

I’m making a design assumption that our quadcoptor will be a square cross configuration and not an, “X-wing” form-factor. I don’t have a particular reason other than wanting to sketch it this way and keep the diagrams and coordinate reference frames as simple as possible.

ref_frames

Angular Equations of Motion

The Bouabdallah paper refers to the Lagrangian method before summarizing the equations of motion for the 3 rotational axes as described above. I had some difficulty following the author’s use of Euler Angles and I’m not sure the derivations use these angles consistently within the paper.

I’m going to try a Newton-Euler method (Morin p.376) because although the Lagrangian method is elegantly simple I find equating the rate-of-change of body momentum to applied forces and torques straightforward. In this case I get less confused between using body angles and Euler angles (body angles relative to the fixed-frame. I hope to repeat the derivations using the Lagrangian at some point.

Recommended Reading

I was referred to the following text as a, “classic”. The first edition was published in 1965, second edition in 1991. I acquired a copy. I see why it is revered. The first 10 pages of the first chapter lay-out the equations of motion for aircraft more concisely than any other text, paper, or online reference I have found. The clarity I need is Blakelock’s labeling and naming of the various body displacements, angles, rates-of-change of position and rotation, linear and rotational momentum, and his description of the Euler angles. In my opinion many other sources are not as clear on these points.

Retired USAF Colonel Blakelock died in 2015. Click his picture to read his obituary. I am thankful he left behind this great book for us.

It is confusing to think back-and-forth from the, “inertial” frame and the “body frame”. Specifically, my mind jumps to the rotations on the body that an eventual gyro package will measure as (p,q,r) where…

p = Rotational rate of the body rolling about the body X-axis.

q = Rotational rate of the body pitching about the body Y-axis.

r = Rotational rate of the body yawing (twisting) about the body Z-axis

Then there are, “Euler Angles” (\Phi , \Theta , \Psi )  representing roll, pitch, and yaw, respectively. These are relative to the fixed-frame: the local NED coordinate system in our case.

In some papers and texts capital and lowercase Greek letters, introduction of other symbols, presence or absence of the (p,q,r) definitions and other information can create confusion and doubt that one is indeed treating the accelerating coordinate frame(s) correctly at the same time accounting for vehicle dynamics.

Blakelock’s first 10 pages offer a clean and clear description that we will later see gets us to the equations in our main reference for the quadcoptor: the Bouabdallah paper.

Blakelock’s straightforward equations-of-motion derivation is included here:

blakelock_eom2

Timeout: Angular Momentum Clarification

In a linear system we know the simple statement, “momentum equals mass times velocity”. We know that linear motion in three dimensions means that velocity is a vector with x-, y-, and z-components. Mass for linear momentum is lumped to a point (our ball-bearing example above). It is a scalar. So, in a linear system we’re left with a momentum vector: the mass times the 3 components of it’s velocity.

A rotational system in three-dimensions will have a rotational momentum vector as well, but in this system, “mass” is represented as, “moment of inertia (I)” and velocity is angular velocity of the entire body represented as a three-coordinate vector \omega or \Omega  (capital or lowercase may vary but it is almost always this Greek letter, “omega”that is used).

At any instant you can think of this \omega-vector as the axle on which the body spins.

To complicate matters further, what was a scalar mass for the linear system becomes not simply a vector of inertia about various axes but a 3×3 matrix!

we’re going to ultimately arrive at a 3-element vector of angular momentum when we multiply our, ‘I’ matrix by our body angular velocity vector \vec\omega .

\vec L\:=\:\b I\vec\omega

First let’s get more comfortable with the inertia matrix, ‘I’.

Platform Inertia Matrix

The Blakelock text above skips over the details of the, “Inertia Tensor”: a matrix that represents to moment of inertia of the body. These are the, “I” terms. Good treatment of this topic is provided by David Morin in Ch. 9 of, “Classical Mechanics with Problems and solutions“. This Chapter is extracted here for your reference. The matrix describing a body’s moment of inertia is covered beginning on P. 376.

MorinCh9

Angular Momentum Changes

Page 393 in Morin above states, “Angular momentum, ‘L’ in the body frame changes relative to the fixed-frame (our NED frame) by two effects”..

  1. Body coordinates change relative to the body frame.
  2. The body frame itself rotates with respect to the fixed-frame.

Morin sets \vec L_0 = \vec L  as the momentum vector of the body at any given instant and imagines painting it on the body at this instant. The rate-of-change of momentum with respect to the fixed-frame is then…

\frac{d\vec L}{dt} = \frac{d(\vec L- \vec L_0)}{dt} + \frac{d\vec L_0}{dt}

That first term is the rate-of-change of momentum with-respect-to the body frame. This is what someone sitting on the rigid body would measure:

\frac{\delta \vec L}{\delta t} = \frac{d(\vec L- \vec L_0)}{dt}

The remaining \frac{dL_0}{dt}  term is the rate-of-change of the body-fixed momentum vector which is the body momentum at this instant as defined above. The momentum vector in this instant is assumed, “painted” on the body like our body axes. This temporarily, “constant” momentum is rotating with respect to the fixed NED frame on the ground at this instant. This is the \frac{dL_0}{dt}  term.

Morin reminds us that the fancy sum above that gives us total change in momentum \frac{d\vec L}{dt} is saying just this:

The total rate-of-change of angular momentum is the rate-of-change of momentum relative to the moving frame PLUS the rate-of-change of the momentum of the moving frame with respect to the fixed frame.

When we have any number of, “reference frames” for vectors this is how we map any vector from one frame to another. There is nothing magic or new here, even though the nomenclature looks fancy.

Body-Fixed Momentum Vector Relative to NED Frame

Just as Blakelock on page 13, Morin p. 93 gives us,

\frac{d\vec L_0}{dt} = \vec\omega \times \vec L

In Morin chapter 9 you can see derivation of the, “Inertia tensor” but it simplifies to a diagonal matrix when we select the body-frame axis to be principal axes of the body.

I =\begin{bmatrix} I_x & 0 & 0\\ 0 & I_y & 0\\ 0 & 0 & I_z \end{bmatrix}

Blakelock doesn’t zero the lower-left and upper-right terms of the inertial matrix. For a fixed-wing aircraft he maintains a cross-coupling of the inertia terms associated with roll and yaw. If you study his figure 1-2 in the book section above you can make sense of this for an aircraft that does not appear symmetrical front-to-back. However, our quadcoptor is symmetrical here too. Our inertia matrix simplifies to the diagonal matrix above.

This leaves us with L_i = I_iw_i for each body axis x, y, and z. The cross-product  \vec\omega \times \vec L  can be computed via the determinant of the matrix below. The determinant of a 3-by-3 matrix can be had by augmenting it with the first two columns as shown and subtracting left-diagonal products from right-diagonal products.

\begin{vmatrix} \hat x & \hat y & \hat z\\ w_x & w_y & w_z\\ I_xw_x & I_yw_y & I_zw_z \end{vmatrix} \begin{matrix} \hat x &\hat y \\ w_x &w_y \\ I_xw_x & I_yw_y \end{matrix}

Below are the hand calculations for the  \vec\omega \times \vec L  cross-product. The black lines represent the, “right diagonal” products. the red lines represent the left-diagonal products.

These hand calculations agree with Blakelock equations 1-32 (the central term with the principal axes moments of inertia in parenthesis). Our rotational rates match Blakelock’s as follows:

\omega_x = P \\ \omega_y = Q \\ \omega_z = R

These are the body rates-of-rotation about the principal axes of the body: the frame, “painted” on the body. We will soon get into details for an Inertial measurement unit (IMU) with gyros that will give us (P,Q,R) data for our platform controller.

Rate-of-change of body momentum relative to the body frame: \frac{\delta \vec L}{\delta t}

This is what we would measure if we were sitting still on the body.

As in Blakelock’s equation 1-32 we will assume the primary-axis inertia terms are constant. For aircraft, missiles, and the like fuel burn changes mass and distribution of this mass within the body. Our battery-powered quad will not experience these changes, but it worthwhile to consider this effect should we desire to deliver a payload with our quadrotor. When we drop the load the mass and moments would change.

Blakelock indicates that even when this is relevant for fuel-thirsty aircraft the above assumption works over time intervals much longer than our tens-of-Hz platform control loop. We would implement a, “sliding” controller or, “gain scheduling” system that effectively remodels our design periodically as fuel is consumed or payload is released. between those slower updates we would take our moment-of-intertia, “I” terms as constants.

Thus, we have

\frac{\delta \vec L}{\delta t} = (\dot PI_x, \: \dot QI_y,\:\dot RI_z)

Combining this with the result above our equations of rotational motion become…

\sum_{}^{} T_x = I_x\dot\omega_x\:+\:(I_z-Iy)\omega_y\omega_z \\\sum_{}^{} T_y = I_x\dot\omega_y\:+\:(I_x-Iz)\omega_x\omega_z \\\sum_{}^{} T_z = I_x\dot\omega_z\:+\:(I_y-Ix)\omega_x\omega_y

If you assume no applied Torques T (natural rotational response of the system) you can rearrange the above equations to match the Bouabdallah paper equation (6). Note that the, “I” subtractions are reversed because assuming zero applied torque we move terms to the other side of the equal sign and solve for the angular acceleration.

These match when we relate Bouabdallah’s terms to our terms above as follows…

\dot \phi = \omega_x \\ \dot \theta = \omega_y \\ \dot \psi = \omega_z \\ \ddot \phi = \dot \omega_x \\ \ddot \theta = \dot \omega_y \\ \ddot \psi = \dot \omega_z \\

Equations of Linear Motion

The Bouabdallah paper describes a test-stand for which various various control strategies are studied. The test stand pins the CG to a point, constraining it to angular motion only. The paper does not cover details for linear motion.

Blakelock chapter 1 above offers good descriptions for modelling linear motion. Below is his conclusion. If you look back at his Equation 1.15 you can see how each primary axis velocity couples into acceleration of the other axes via another cross-product. I’m going to leave this here for now because depending on how we measure position and velocity relative to the body and/or fixed NED frame these terms may or may not be relevant in our design.

 

Concluding Remarks…

It looks like we have a good handle on the equations of motion now. I am more confident in nomenclature and with the mental gymnastics involved with thinking within the two reference frames:

  • The North-East-Down (NED) frame pinned to the earth at our take-off point.
  • The Body frame we can think of as permanently painted along primary (x, y, z) axes of our quadcopter.

I’m going to pause here, and pick-up the applied torque and force elements of the equations-of-motion in the next post. These are the following terms in the equations above.

\sum_{}^{} T_{(x,y,z)}\\ \sum_{}^{} F_{(x,y,z)}\:

They are the torques and forces that will act on the quadcoptor. They will include the desired outputs from the propulsion units, and we can also model external factors such as wind. We will need a controller design that is robust to some applied, “disturbances” from wind.

I’m looking forward to getting into sensor and control system designs!




Quadrotor Dynamic Model: Propeller Gyroscopic Effect.

I copied the equations below from the last post, where we see we matched our Bouabdallah paper‘s equations for rotational motion. We know some propeller input is going to be a part of the applied torque. There is going to be a set of equations for linear motion too, but let’s clarify what is going on here regarding torques about our 3 body axes. Let’s look at these equations and…

Review Rigid body rotation gyroscopic effect

\sum_{}^{} T_x = I_x\dot\omega_x\:+\:(I_z-Iy)\omega_y\omega_z \\\sum_{}^{} T_y = I_x\dot\omega_y\:+\:(I_x-Iz)\omega_x\omega_z \\\sum_{}^{} T_z = I_x\dot\omega_z\:+\:(I_y-Ix)\omega_x\omega_y

Remember Newton’s Second Law \sum F_i=m\cdot a_i? This tells us that the sum of the forces on a mass in the i’th direction must equal the mass of the body times the acceleration in that direction.

If we take one of the equations above and rewrite it we can recognize this familiar form, albeit in rotational nomenclature. The rightmost term is the gyroscopic effect of the entire rigid body. It tells us that if there is rotational rate about y- and z-axes a torque is induced about the X-axis. If we neglect this term you see we’d be left with  \sum_{}^{} T_x = I_x\dot\omega_x in simple second-law form.

It’s a little confusing, but all we’ve done is already offer one item to the, “sum-of-torques”: the gyroscopic effect resulting from the rigid body rotation. Notice when we moved the term in parenthesis to the other side it looks like an added torque now but we swapped the subtraction in parenthesis so the math is the same.

I_x\dot\omega_x = \sum_{}^{} T_x \:+\:(I_y-I_z)\omega_y\omega_z

Gyroscopic effect of Propellers rotation.

“tail-dragger” airplane example

I made a video below to describe a phenomenon familiar to pilots, particularly pilots of powerful, “tail draggers” like the P-51 mustang(below). A propeller is a large spinning disk. when we, “tip” the plane of rotation of the propeller: pitch the airplane, a torque on the body results due to this gyroscopic effect.

Imagine you are in the cockpit of the P-51, speeding down the runway (lucky you). That monster prop has a healthy angular momentum vector pointing perpendicular to the prop. Your tail is dragging until you get enough lift from your wings (due to streamline curvature, if you recall a recent post!). At some point your tail will lift quickly and you’ll be parallel to the ground before your front wheels leave the ground. This looks like a quick, “pitch” of say 15-degrees based on the photo above. Let’s say it took a second for the tail to rise as we’re nearing take-off, so the pitch rate is 15-degrees per-second as the tail wheel leaves the ground.

That same pitch is pitching the, “plane” of the spinning propeller forward, or more to the vertical. When you watch the video below and inspect the equations below you will see that this, “pitch forward” of the propeller at the above angular rate of the airplane going from dragging-to-level results in a, “gyroscopic torque” about the airplane’s vertical axis: you feel the plane want to, “yaw” or twist so you’ll need to correct with some rudder.

The effect can be negligible

Our propellers are horizontal in the body-plane, but there is a torque induced on the body that results from propeller rotation and rotation of the propeller plane about the x- and y-axis of our system too…at least in theory. It may be negligible depending on propeller moment of inertia relative to other effects, but let’s model it anyway.

We can choose to include it or neglect it later when we perform a sensitivity analysis. This is where we can look at rough magnitude of various effects and decide to include them in our model or cancel the terms in the equations: ignore them as too small to make a difference. We’ll do this systematically and not guess at this point. Besides, this is a neat phenomenon for us to appreciate in detail.

Video Explaining Equation Derivation

Gyroscopic effect of our propellers

The following PDF derivations arrive at the same terms as the Bouabdallah paper (equation 8). There’s a sign convention difference between the paper and my derivations, but the result is the same:

T_{x,prop}\:=\:I_{prop}\cdot\dot{\phi}(\omega_2+\omega_4-\omega_1-\omega_3)\\ T_{y,prop}\:=\:I_{prop}\cdot\dot{\theta}(\omega_1+\omega_3-\omega_2-\omega_4)

The first equation tells us that there will be a torque about our body X-axis when there is angular rate-of-change to our pitch: \dot \phi . This means our quadcopter will feel a  body roll-inducing torque when it pitches. This is due to the gyroscopic effects of the 4 propellers.

The second equation tells us that there will be a torque about our body Y-axis when there is angular rate-of-change to our roll: \dot \theta. This means our quadcopter will feel a  body pitch-inducing torque when it rolls.  Again, this is due to the gyroscopic effects of the 4 propellers.

The hand calculations below arrive at the above equations following the process in the video. The video above explains how to go about the hand calculations for a single propeller.

Hand Calculation PDF Notes

RB 2018-07-11 17.27.01

 

Conclusion

Above we have one more term for our sum-of-torques about the body X- and Y-axes. The plus-minus signs represent our set-up with two clockwise spinning propellers and two counter-clockwise spinning propellers. If you expand these equations you’ll see I just performed the steps in the above video Eight times and combined the equations. We need eight steps because we have 4 propellers and we need to account for the plane-of-rotation rate-of-change about both the x- and y-axis. Rotation about the body z-axis contributes no gyroscopic effect from the propellers because it is just like twisting the axis of the propeller. The plane-of-rotation of the propeller does not change at any rate in this case so no torque is induced on the body.

T_{x,prop}\:=\:I_{prop}\cdot\dot{\phi}(\omega_2+\omega_4-\omega_1-\omega_3)\\ T_{y,prop}\:=\:I_{prop}\cdot\dot{\theta}(\omega_1+\omega_3-\omega_2-\omega_4)

We’ve now accounted for gyroscopic effect of the rigid body rotations and the propeller rotations.

Next we’ll add the propeller thrust components. These will be our inputs to counter-act the above effects and unknown disturbances (wind) in order to maintain the attitude of our quadcoptor and move it where we want it to go.




Finalizing Equations Of Motion: Thrust Inputs from Propellers

This post explains how we determine propeller thrust and drag factors for our quadcopter project.

The last couple posts have been working-out the sum-of-torques on our quad-copter. The first, “unforced” model considers the gyroscopic effect of the total air-frame. That was the first, “equations of motion” post a few weeks ago.

Next we looked at the torque induced by the gyroscopic effects of the spinning propellers in the last post.

Now we’ll consider the propeller drive inputs.

Propeller Drive Control Input

I started the, “equations of motion” with the sketch below, so we’ll use it again here as we model the command inputs from the propellers: F_{1-4}.

ref_frames

Torque About Body X-axis

Assume the length of the arm from the center of the quadcopter to the propeller is, ‘l’ then…

\tau_x=F_2\cdot l - F_4\cdot l

The, ‘F’ terms are thrust forces. We know from an earlier post that this thrust from a propeller is going to be a function of propeller speed squared: \Omega^2:

F = b\cdot\Omega ^2

Where, ‘b’ is a, “thrust factor” we need to determine. After you watch the video below you’ll understand how we get…

b =\frac{C_t\cdot\rho\cdot D^4}{4\pi^2}

where…

C_t = experimentally \ determined \ propeller \ thrust \ coef.\\ \rho \ \ = \ air\ density.\\ D=propeller\ diameter.

our torque about the X-axis equation above becomes…

\tau_x= bl(\Omega _2^2-\Omega _4^2)

Torque About Body Y-axis

\tau_y=F_1\cdot l - F_3\cdot l

By the same logic above for the X-axis torque, this Y-axis torque becomes…

\tau_y= bl(\Omega _1^2-\Omega _3^2)

Torque About Body Z-axis

The torque about the Z-Axis is proportional to the squared angular velocities of the propellers. We subtract the squared counter-clockwise props from the clockwise props to lump them. The proportionality constant (the multiplier) is going to be a, “drag factor”: d.

See the video for an explanation on how we get ‘d’ from the experimentally determined, “power coefficient” for a candidate propeller.

\tau_z = d\cdot (\Omega _2^2+\Omega _4^2-\Omega _3^2-\Omega _1^2)

d =\frac{C_p\cdot\rho\cdot D^5}{8\pi^3}

C_p = experimentally \ determined \ propeller \ power \ coef.\\ \rho \ \ = \ air\ density.\\ D=propeller\ diameter.

Deriving the, ‘b’ and, ‘d’ terms above

‘b’ is our thrust factor and, ‘d’ is our drag factor. They are arrived at through experimentally determined thrust and power coefficients for our propeller.

An excellent resource is offered by the University of Illinois and Urbana-Champaign: the UIUC Propeller Data Site. There are wind-tunnel data for a variety of hobby-craft propellers. The following video explains how we get our, “factors” for our equations of motion from the experimental, “coefficients”…

Propeller Thrust Coefficient

I am going to estimate my thrust coefficient from the, “static” data because a quadcoptor is not going to travel very fast along the axis of the propellers. I will likely select a propeller model included in the UIUC data sets. My preliminary estimate for the static thrust coefficient after looking at a number of experimental data sets is…

C_t =0.1

From this I can calculate, ‘b’ as described above.

Propeller Drag Coefficient

Remember the motor+gearbox+propeller post from a few months ago (already!)? There I was looking at now to map the torque load of the propeller back through the gearbox to the motor shaft to get all the parameters correct in modelling the DC motor drive. In that post I write…

“Aerodynamic torque for rotor R:    D = 2.91\times10^{-7} \frac{Nms^2}{rad^2}    (page 9, Gafvert) for what looks like a comparable hobby-sized propeller lab set-up”. Until now, I had no idea how this term was derived.

Confusion…

It was a challenge for me to find a propeller, “drag coefficient” explicitly stated in my references. Then I read, “While an airfoil can be characterized by relations between angle of attacklift coefficient and drag coefficient, a propeller can be described in terms of advance ratiothrust coefficient, and power coefficient” here. That statement cleared some fog.

Next I found an online question from a person facing the same confusion I was experiencing. The question married the above quote to the same UIUC data I was looking at. The dialog there ultimately gets to a torque term resolved from the UIUC, “power coefficient”.

From there I made sense of it and as expressed in the video above. I derive my, ‘d’ term from the experimental, “power coefficient” as explained in the video. As the video explains, I’m estimating…

C_p=0.7

For the APC Carbon Fiber 7.8-inch propeller, if we plug-in numbers we calculate, ‘d’ as…

d =\frac{C_p\cdot\rho\cdot D^5}{8\pi^3}=\frac{(0.7)(1.23\frac{kg}{m^3})(0.175m)^5}{8\pi^3}=5.65\times10^{-7} Nm s^2

This is close to the number from page 9 of Gafvert above. I don’t know what specific propeller he is referencing, but we know it is in the family of these hobby propellers. I’m satisfied with this method for arriving at our drag factor from the UIUC data.

Conclusion…

We have our thrust and drag, “factors” for our propellers (in control-system parlance we can think of these as, “gain” terms).

When we get into the design iterations we’ll likely try other propellers with different coefficients. We can revisit the UIUC data and compute a range of, ‘b’ and, ‘d’ terms as we did above. This is how we can assess propeller selection impact on performance of our quadrotor.