Dynamics of the centrifugal governor : Lagrange method with Vex simulation example

June 18, 2015

This post derives the equations of motion for a motor driving a variable centrifugal load similar to what one would see on the classical steam engine centrifugal speed governor first invented by James Watt. (Fig 1) Source(https://upload.wikimedia.org/wikipedia/commons/c/c7/Boulton_and_Watt_centrifugal_governor-MJ.jpg)

Non-Feedback Dynamics

We are only concerned here with the response of the engine/centrifugal load without any feedback to the engine from the governor.   The analytical diagram for the governor without feedback is:

centrifugal governor analysis diagram3

Typically a motor with a large inertial load has a first order exponential  response with a time constant

tau = w_free*I_load/torq_max*gr^2

where w_free is the motor free speed , torq_max is the stall torque of the motor and gr = gear ratio of output shaft speed / motor speed.

The centrifugal weights  change the moment of inertia I_load with output shaft speed hence giving the motor a variable time constant.    One can bound the motor response between two exponential responses… one with a tau_min and the other with tau_max where min and max tau are respectively associated with the min and max inertial load which are generally set by the mechanical limits of the governor.   The motor current response is nonlinear and  more gradual than if operating with the full inertia load.   This in turn can keep the motor protection PTC fuses from tripping without explicitly controlling the current through the motor controller.

I will post a simulated response later for a Vex system but first we need the equations of motion.

To find the equations of motion I decided for my review to use the Method of Lagrange.  For a given system one needs to compute the Lagrangian, L, from the kinetic energy, T and the potential energy V.

L = T – V

We then obtain a system of differential equations from the following computation for each generalized variable q_i.

d(∂L/∂q_i_d)/dt = ∂L/∂q_i + Q_i

where q_i_d  = d(q_i)/dt  and Q_i are the generalized non conservative forces doing work on the system along the q_i direction.     This is a cook book formula that is often times simpler to use than Newtons force equation F = ma.     The resulting equations will be the same with either method.   Many good references on Lagrange method  can be found with a simple Google search.[reference Wikipedia Euler Lagrange Method]

The  generalized variables of interest are :

q1 = psi , the angle that the motor is turning through.

q1_d = w_axel  = dpsi/dt , the speed at which the motor is turning the governor axles.

q2 = alpha, the angle of the arm from vertical.

q2_d = w_arm = dalpha/dt , the speed that the arms are moving up and down relative to the vertical.

The mass of the weight on one arm is m and the length of the arm is l.   Moments of inertia relative the pivots are I_axel and I_arm.    I_arm is constant and I_axle varies due to the movement of the arms.

T = 1/2 * I_axel * (w_axel)^2  +  1/2 * I_arm * (w_arm)^2

V = m*g*l*(1-cos(alpha))

I_axel = m*(l*sin(alpha))^2

I_arm = m*l^2

Substituting for the generalized variables

L = T-V = 1/2*m*(l*sin(q2))^2 *q1_d^2 + 1/2*I_arm*(q2_d)^2 – m*g*l*(1-cos(q2))

∂L/∂q1 = 0;

∂L/∂q1_d = m*l^2*sin(q2)^2*q1_d

∂L/∂q2 =  m*l^2*sin(q2)*cos(q2)*q1_d^2  – m*g*l*sin(q2)

∂L/∂q2_d= I_arm*q2_d

Now compute d( ∂L/∂q_d)/dt and set equal to dL/dq + Q  where Q are the generalized forces that are non conservative associated with each q.

For q1:

d(m*l^2*sin(q2)^2*q1_d)/dt = m*l^2*(2*sin(q2)*cos(q2)*q2_d*q1_dot + sin(q2)^2*q1_dd)  = torq

1) q1_dd = -2*cos(q2)*q2_dot*q1_d/sin(q2) + torq/(m*l^2*sin(q2)^2)

For q2:

I_arm*q2_dd =  m*l^2*sin(q2)*cos(q2)*q1_d^2  – m*g*l*sin(q2)

or

2) q2_dd = ( cos(q2)*q1_d^2 – g/l)*sin(q2)

Steady State requirements

In the steady state, q1_dd = q2_dd = q2_dot = o.   This then gives the  q2 as a function of q1_d or alfa as a function of the axle speed , w_axle.

From equation 2    cos(alfa) = g/(l*w_axle^2)

Since cos(alfa) < =1     the minimum axle speed  w_axle_min >=sqrt(g/l)    Notice that this requirement is independent of the centrifugal mass very much like a pendulum period is independent of mass and only a function of the pendulum length , l .  w_axle_min is  the minimum speed before the arms move.

Vex Example

In the Vex centrifugal experimental device shown here , the arms are about 2.5 inches in length so we would expect w_axle_min = sqrt(32.3*12/2.5)= 12.4 rad/s or 118.5 rpm.   There appears to be a gearing between the motor and axle of 7 : 1 so the minimum speed of the motor is  118.8/7 =17 rpm   which means that the motor driving  can be a Vex 393 with standard gearing (100 rpm max).      However, there is a mechanical minimum of around 45 degs   which means the minimum speed to move is   17 rpm/sqrt(cos(45 deg)) = 20.2 rpm.   The speed to fully extend the  arms (alpha = 75 deg ?) would be  17 rpm/sqrt(cos(75 deg)) = 33.4 rpm.

System Time constant with alfa = 45 deg and 75 deg. 

tau = w_free*I_load*gr^2/torq_max  

I am guessing that one arm has a mass of about 1 oz  0r .0283 kg

I_45 = m*(l*sin(alfa))^2 = .o25*(2.5in*.707*.0254 m/in)^2 = 5.038 E-05 kg m^2

I_75 = I_45*(sin(75)/sin(45))^2 = I_45*1.866

tau_45 = w_free*I_45*2*gr^2/torq_max   (note: factor of 2 added for two arms)

Since torq_max = 1.67 nm, w_free = 100 rpm*6.28/60 = 10.5 rad/sec , gr = 7

tau_45 = 10.5 * 5.038E-05*2*49/1.67 = .037sec

tau_75 = tau_45*I_75/I/4=  .069 sec 

Simulated response

The differential equations were implemented in a simulation to show the time response to a step speed input.  The data was generated by my excel program.   The red line in the time response shown below is the arm angle starting from 45 deg and quickly moving to 75 deg as the speed increases.    This expansion occurs so quickly that most of the motor response is with the maximum moment of inertia.   In the second figure, the %current error is shown with fixed and varying moments of inertia.

Step Response Centrifugal governor

As you can see after about .o2 seconds, the motion of the governor is limited at its maximum alfa = 75 deg and the variable inertia current response (blue curve) is close to the slower fixed inertia response (yellow curve) for most of the current response.

%current response with centrifugal governor1The o

In this particular experiment, the objective was to keep the current response below the response with the maximum inertia to keep the stress of the motor.  In this case the objective was not met.    The mass movement was too fast with speed so the governor limited shortly after the initial speed increase.    Increasing the length of the arms expands the time scale but doesn’t change the early saturation.   One solution would be to add spring tension to the expanding weights so that a higher speed would be required to saturate.    This will be addressed later by adding a spring potential energy to the Lagrangian and deriving a new set of equations.


Vex Note: How a single flywheel ball shooter minimizes the effect of ball mass variations

May 28, 2015

Nothing but Net 2015/2016 competition game involves shooting 4 inch balls that can have a 10% variation in mass.    We know that trajectory range ,R = V^2/g*sin(2*theta) so it  is dependent upon the square of the ball release speed , V, and shooter elevation, theta.   Mass does not enter into the equation unless it affects V.

Ball release energy :

Suppose we use a Vex 5″ diameter wheel as a flywheel and rotate it a w_wheel angular speed.      As the ball leaves the shooter, it will have a V = r_wheel*w_wheel/2.  e.g. half of the flywheel tangential speed.    The ball will have a spin rate , w_ball = V/r_ball.    The energy of the ball, E_b , is the sum of the ball translational energy and rotational energy.

E_b = 1/2*m_ball*V^2 + 1/2*I_ball*w_ball^2

where I_ball = 2/5*m_ball*r_ball^2 (solid sphere of uniform density).

so Eb =  1/2*m_ball*V^2( 1+2/5)  .   (corrected 5/29 Was 1/2*m_ball*V^2( 1+4/5)  So the rotational energy adds  40% more to the translational energy.  Rewriting in terms of w_ball gives

E_b = .7*m_ball*w_ball^2*r_ball^2  

Wheel Energy:

E_wheel = .5*I_wheel*w_wheel^2.  where

I_wheel = m_wheel*(r_wheel*.84)^2  (ref blog post https://vamfun.wordpress.com/2015/05/17/finding-the-moment-of-inertia-of-a-vex-wheel-using-parallel-axis-theorem/)

Energy Conservation:

E_wheel_initial = E_wheel_final + E_ball     This assumes that the wheel is not being powered by the motor during launch and that the extra energy needed for the ball comes from the flywheel.   Also, friction and ball compression energy losses are assumed zero to simplify this analysis but can be significant in actual percentages derived.   I am focusing  on how increasing flywheel mass lowers the percentage range errors caused by ball mass variations.

E_wheel_initial/E_wheel_final = (1 + E_ball/E_wheel_final)

Lets expand E_ball/E_wheel_final

E_ball/E_wheel_final = (.7*m_ball*w_ball^2*r_ball^2)/(.5*I_wheel*w_wheel_final^2)

= 1.4*m_ball*w_ball^2*r_ball^2/(m_wheel*r_wheel^2*.84^2*(2*w_ball*r_ball/r_wheel)^2)

= .4954*m_ball/m_wheel

SInce  m_ball = 60 g and m_wheel = 180 g   m

_ball/m_wheel = 1/3

So  E_ball/E_wheel_final = .165    for a single 5″  wheel flywheel     .165/n for n flywheels.    So the ball energy is almost equal to the 1/6 final energy of the wheel

Range Tolerance analysis:

So how does R vary with m_ball from all this.   Well , we know the range is proportional to V^2 which is proportional to w_wheel_final^2 which is proportional to E_wheel_final.

From above E_wheel_final = E_wheel_initial/(1+ .4954*m_ball/m_wheel)

So due to proportionality of R and E_wheel_final we can say

R/R_0 = ((1+ .4954*m_ball_0/m_wheel)/(1+ .4954*m_ball/m_wheel))

where R_0 and m_ball_0 are the nominal values without errors.

We can use R range= R_0(1+ %e_r)   and m_ball = m_ball_0*(1 + %e_m_ball) to work with % changes.

Then with some manipulation we can get %e_r as a function of %e_m_ball

%e_r  = -%e_m_ball/(2.02*m_wheel/m_ball_0 +1 + %e_m_ball)

Now m_wheel = n*.180 kg   and m_ball= .06 kg  so we can write an approx.

%e_r = -%e_m_ball /( n*6.06 +1)     where n is the number of 5″ vex wheels.

Lets put in a few numbers:

Assume %e_m_ball = 10%  then the range error is

n = 1, %e_r =  -1.42%

n = 2, %e_r =  -.76%

n = 3, %e_r =  -.52%

n = 4, %e_r =  -.40%

n = 5, %e_r =  -.32%

So you see the benefits of having a higher  flywheel mass to ball mass ratio.   The use of  two 5″ wheels in a single wheel design can reduce a potential 10% range error from ball mass variations  to  1% ( less than a ball radius).   To keep the spin up time to a reasonable number of seconds requires about 2 393 motors per wheel so 2 wheels costs 4 motors.   So there is a motor tradeoff to get that  higher accuracy with heavier flywheels.


Finding the moment of inertia of a Vex wheel using parallel axis theorem

May 17, 2015

vex 5 in wheel

The new Vex game “nothing but net” might involve rotating shooter wheels.  We know that if all the mass of a wheel was located on its rim then the moment of inertia about its rotating axis  (I_rim) would be

I_rim = r^2 * m   where m is the mass of the wheel and r = radius of wheel.   But we know that the wheels actually  have mass that is unevenly distributed along the radius so the moment of inertia I_wheel will be less than I_rim.

Easy experiment to determine I_wheel if we know its mass.

We can determine I_wheel experimentally using the parallel axis theorem and the dynamics of a pendulum.

Parallel axis theorem says that any object that is rotated about an axis parallel to and a distance , d, from an axis going through the centroid of the object will add an amount =  m*d^2 to the moment of inertia about its centroid.  I.e.

I_parallel = I_centroid + m*d^2 .

Suppose we now swing the mass, m,  about the parallel axis like a pendulum  using just the torque from gravity pulling on the mass.      It is easy to show that the period, T , of the pendulum is related to the distance , d, and the moment of inertial , I_parallel, by the following formula:

T = 2*pi*sqrt(I_parallel/(d*m*g))     .   g is gravitational constant and assumes swing angles smaller than say 10 degs from the lowest point on the pendulum path. 

If we measure the period of the pendulum we can rearrange the equation and find I_parallel

I_parallel = T^2*d*m*g/(2*pi)^2 

Once we have I_parallel, we can now use the parallel axis theorem to determine I_wheel.

I_wheel = I_parallel –  d^2*m  =  d^2*m * (T^2*g/d/(2*pi)^2 -1)

(This assumes  that the string has negligible mass relative to the mass of the wheel)

Vex 5 in Wheel experiment:

Given….r_wheel = 2.5 inches ,

wheel mass,   m = 180 gm (0.180  kg)

The pendulum is created by suspending the wheel with a thread 2.75 inches from its center so

d =. 07 m (approx. 2.75 inches)

The average period  T = .668 sec

I_wheel  = d^2*m*(T^2/d*9.8/(6.28)^2 -1)

     = d^2*m*( .248*T^2/d -1)

     = .07^2*.180*(.248*.668^2/.07 -1) = 0.00051 kg m^2

Equivalent radius with a rim only mass r_e

r_e = sqrt( I_wheel/m) = .0533 m ( 2.1 inches)  

This  means that the wheel behaves as if the mass if located at 84 %  of the radius of the wheel which one could almost guess by looking at it.


60 in double reverse 4 bar linear linkage example

June 8, 2014

Here are some pictures of how a double reverse 4 bar using Vex 35 hole link arms. The main challenge is constructing a stable geometry. The lift has a middle link gear train to move the upper 4 bar in sync with the lower 4 bar.

Torque

Approx torque (inlb) required = 2*l_arm*(W_payload_lbs + W_lift_lbs + W_manipulator_lbs)*cos(angle)

height_delta = 2*l_arm*( sin(angle_finish) –  sin(angle_start))

where l_arm is the distance between pivots of the arms in inches, W_payload_lbs is game piece weight, W_lift_lbs  is the weight of the 3 arms used in the lift , W_manipulator_lbs  is the weight of the gripper attached to the upper arm and angle is defined relative to the horizontal.  The max torque occurs when the arms are horizontal or angle • 0 .

If you want to lift 4 1 lb Skyrise cubes with a 3lb lift weight , a 1 lb  gripper and a l_arm of 16 in you would need about 256 in lbs of torque.

As a rule of thumb I use 6 inlb of torque per motor for sizing  the number of motors.  This assumes 2  inlbs of elastic support , 1 inlb of friction and 3 inlbs from active current (about .9 amps) hold per motor used.   With these assumptions a 10:1 gearing and 4 393 motors might do the job.    The lift would take about 5 seconds to go full travel.  I’ll show a more exact torque derivation later.

Height change

With an  angle_finish = 75 deg and the angle_start= -60 deg  a 16 in arm reaches 59 in.  A 16.5 in arm reaches 61 in.   With a chain bar you could  clear an in or two more.20140608-214047.jpgSo 20140608-214101.jpg 20140608-214114.jpg


Getting my feet wet with Arduino Uno

July 18, 2013

Arduino Uno

I have tried to stay away from Arduino mainly because I didn’t want to learn a whole new culture. However, my AlgalitaROV project uses the OPENROV electronic architecture which employs a Cape that uses the Arduino Uno processor which in turn attaches to a BeagleBone Black(rats…another processor culture to learn).
Downloaded and installed the Arduino SDK 1.0.5 from arduino.cc and read through the getting started and reviewed the types of commands it uses. I bought a Arduino Uno starter kit from amazon ($32) which has a breadboard and some jumper wires. plugged in the Arduino Uno and it installed automatically without updating the drivers as suggested in the instructions. First observations:

1) the general purpose I/O (GPIO) are 5 volt at the breakouts. This is great since it is compatible with Vex parts. Current is o limited to 40ma per pin.

2) it can be powered from the USB cable or an external supply that must range from 6v to 20v with a recommendation of 7v to 12v.

3) Software wise in the Arduino you write a void setup() function and a void loop() function. These functions are called automatically by the operating system. setup() runs once and the loop() runs continuously…meaning the hidden main() function contains a while(true){loop()} structure that runs after the setup() runs. By contrast the Vex RobotC template calls a function pre_autonomous() once and executes separate tasks autonomous() and usercontrol(). usercontrol contains an explicit while(true){} which the programmer can change. Not sure that Arduino gives us that capability.

3)There are plenty of software examples and the language is close to C. The only function that seemed misnamed was analogWrite(pin,duty). This doesn’t write an analog output to a d/a hardware as one might assume but puts a 5 volt PWM signal output on a PWM capable GPIO pins. On the Uno 6 of the 14 GPIO pins are PWM capable (3,5,6,9,10,11). Not sure why they didn’t just call it PWMWrite().

4)Servo library is available for use but be careful with servo.write(int val).

Here is the function code from the servo library.  Note that if val is >544 then the command is assumed to be in microseconds rather than servo degrees.  Also note if it is less than 544 that the degree values are truncated to valid servo range of [0,180].    So a val = 255 will produce a 180 degree command.

void Servo::write(int value)
{
  if(value < MIN_PULSE_WIDTH)
  {  // treat values less than 544 as angles in degrees (valid values in microseconds are handled as microseconds)
    if(value < 0) value = 0;
    if(value > 180) value = 180;
    value = map(value, 0, 180, SERVO_MIN(),  SERVO_MAX());
  }
  this->writeMicroseconds(value);
}

First Script.
We need an r/c brushless motor driver. The ESC needs a standard servo drive signal. So I hooked up a Vex pot to the board and read its value [0, 1023] and scaled it to [0,180] which is given to a canned servo program object which has a servo.write which generates the proper servo output to a specified pin.  A motor is a continuous rotation servo  so 0 value represents max reverse speed and 180 value gives max forward speed.   I also included a serial printout of the value to the laptop console.
Hardware Setup:
[click on picture to expand to full size]

Shows breadboard wiring , motor 7.2v battery and external board 5 volt power supply

Shows breadboard wiring , motor 7.2v battery and external board 5 volt power supply

It is important to note that the motor is powered by a separate battery since it draws more current than the board can supply.

Closeup of arduino wiring

Closeup of arduino wiring

Full test setup

Full test setup

Mymotorcontrol.ino

// Controlling a servo/motor position using a potentiometer (variable resistor) 
// by Chris Siegert:  vamfun@yahoo.com   https://vamfun.wordpress.com
// adapted from  Michal Rinott servo code  
// connect the center pin of a potentiometer (pot) to the potpin.  Connect the other two pins to the arduino 5v and ground pins.
// connect the white wire of the motor pwm cable to the servopin.  
// Important!! Connect the red and black wires to an external motor power source. 
// The arduino is not designed to deliver enough current to drive a servo/motor.
#include <Servo.h> 

Servo myservo;  // create servo object to control a servo or motor. 
int servopin = 9; //connects servo to pin 9 which is a pwm pin.
int potpin = 0;  // analog pin A0 used to connect the potentiometer
int val;    // variable to read the value from the analog pin 
int cmd = 0; // cmd to motor [0 to 180]
int cnt = 0; // program cycle counter 
int cnt_max = 33;   // this is the number of delay cycles between serial printouts. (15ms per cnt) 33 counts gives about 2 per second)  
void setup() 
{ 
  Serial.begin(9600);   //start up the serial data port at 9600 baud 
  myservo.attach(servopin);  // attaches the servo on servopin number to the servo object void loop() 
{ 
  val = analogRead(potpin);            // reads the value of the potentiometer (value between 0 and 1023) 
  cmd = map(val, 0, 1023, 0, 180);     // scale it in degrees of servo rotation(value between 0 and 180) 

  myservo.write(cmd);                  // sets the servo position according to the scaled value in deg 
  cnt++; 
  if(cnt > cnt_max)                  // this sets time between serial print outs to cnt_max cycles  or about (15 ms*cnt_max)
  {
    Serial.println(cmd);  
    cnt = 0;      // debug value
  }
  delay(15);                           // waits for the servo to get there 
}

My kind of kitty kat…. no litter required!!

June 24, 2013

“Thanks to its legs, whose design faithfully reproduces feline morphology, EPFL’s 4-legged ‘cheetah-cub robot’ has the same advantages as its model: It is small, light and fast. ”  Quoted from phys.org newsletter.

catrobot


Note: Using a Line Sensor to Update Gyro Heading

February 28, 2013

Using Line Sensor to Update Gyro Heading

The programming challenge for the 2013 Sack Attack game was attempted by team 1508 LancerBots fully autonomous robot. It uses gyro information to generate a navigation solution and then automatically sequences through several way points to solve the route sequences. Unfortunately as mentioned in earlier posts, the gyro drift causes heading corruption which then causes sacks to by slightly missed. One possible solution I am proposing is to use line sensors to update the gyro.
Figure 1 shows the geometry and formulas required to accomplish the update for a special case when the line is parallel to the truth x axis.    The coordinates used by 1508 were specifically chosen so all lines were parallel to either the x or y  axis.  This simplifies the update equations.

Simply put… when a line sensor event occurs, note the position error  with respect to the line (using the Nav coordinates) and divide it by the range from the origin.    This ratio is the gyro drift error in radians.

Error Contributions:  

Errors in the update equation can stem from errors in  the Nav initial condition (<.25 in, 1 deg heading) , encoder scaling (typically <1%) , encoder resolution(< .03 in) , line location measurements (<.25 in), line tape width (<.5 in) and Navigation solution numerical errors (sampling, round off).

The dominant error is the encoder scaling which can add up when the encoder has moved several hundred inches.   Lines can also be use to update the Nav  position solution so that these  errors remain less than 1% distance  from the field origin and not the total distance traveled by an encoder.   This error translates to about .6 deg of heading error from a gyro update.    So we should not expect much better than that from our updates.

General Case: Lines skew to the coordinate axes

To complete this note,  lets consider a more general line that is described by

c = a*x + b*y

where a,b and c are constants and again x and y are truth coordinates.   Now we have three unknowns (x,y and delta) and three equations:

x’ = x*cos(delta) + y*sin(delta)

y’ = -x*sin(delta) + y*cos(delta)

c = a*x + b*y

First we solve for x and y from the first two equations which is a simple inverse:

x = x’*cos(delta) – y’*sin(delta)

y = x’*sin(delta) + y’*cos(delta)

Make the small angle substitution gives

x = x’ – y’*delta

y = x’*delta + y’

Now substitute x and y into the line equation

c = a*(x’ – y’*delta) + b*(x’*delta + y’)

Now solve for delta

c -a*x’ – b*y’ = (b*x’ – a*y’)* delta

or

answer:

delta = (c -a*x’ – b*y’ )/(b*x’ – a*y’)

Check the special case above :  y = c, a = 0, b = 1

delta = (y – y’)/(x’)  : checks ok