Launch 11 Wind Tunnel

We treated Launch-11 as a wind tunnel test by sweeping the fins back and forth as the rocket flew through the air. Hopefully we can get some useful information from it.

We had two principle sensors on the rocket. An accelerometer and a rate gryo. We can differentiate the rate gryo to get angular acceleration, and integrate the accelerometer to get velocity and height. This along with fin position should give us enough information to solve for $C_L$

Torque

Fundamentally roll acceleration should be based on torque and mass moment of inertia. We have a good guess for the rocket moments, so we believe we can measure torque.

$$\begin{equation} \tau = I\ddot\theta = \mathbf{r} \times \mathbf{L} \end{equation}$$

Where

$\tau$ is the torque on the rocket
$I$ is the mass moment of inertia in the torque axis
$\ddot\theta$ is anguar acceleration
$\mathbf{r}$ is the vector from the center of the rocket spin axis to the center of lift of a wing
$\mathbf{L}$ is lift from a wing

Now we can make the simplifying assumtion that lift is perpendicular to $\mathbf{r}$ so we ingore the whole vector thing and just mulitply. And we can expand our lift equation (note: we have 4 canards!):

$$\begin{equation} I\ddot\theta = 2C_{L\alpha}\rho v^2 A r \end{equation}$$

Solving for $C_{L\alpha}$:

$$\begin{equation} C_{L\alpha} = \frac{I\ddot\theta}{2\rho v^2 A r} \end{equation}$$

Do we have enough information from our two sensors?

On the right hand side we have

$I$: measured ahead of time, geometry and mass dependent.
$\ddot\theta$: differentiate roll gyro
$\rho$: double integrate accelerometer, then use atmosphere model
$v$: integrate accelerometer
$A$: measured ahead of time, flight independent
$r$: measured ahead of time, flight independent

Since we recoreded $\alpha$ directly we should have everything we need.

Let's load the Launch 11 data:



In [1]:

    
%matplotlib inline
import numpy as np
import utils
from matplotlib import rcParams
rcParams.update({'font.size': 14})


columns = np.loadtxt("data/ADIS.csv", delimiter=',', unpack=True)

seqn = columns[0]
timestamp = columns[1]
gyr_x = columns[3]
gyr_y = columns[4]
gyr_z = columns[5]
acc_x = columns[6]
acc_y = columns[7]
acc_z = columns[8]
mag_x = columns[9]
mag_y = columns[10]
mag_z = columns[11]

# normalize time to liftoff, put in seconds
t_0 = 202945201730391 # nanoseconds in FC time frame
timestamp = np.subtract(timestamp, t_0)
timestamp = np.divide(timestamp, 1e9)

# limit to just boost and coast data (don't care after the chutes open
timestamp = timestamp[1700:30000]
gyr_x = gyr_x[1700:30000]
gyr_y = gyr_y[1700:30000]
gyr_z = gyr_z[1700:30000]
acc_x = acc_x[1700:30000]
acc_y = acc_y[1700:30000]
acc_z = acc_z[1700:30000]
mag_x = mag_x[1700:30000]
mag_y = mag_y[1700:30000]
mag_z = mag_z[1700:30000]

Roll Rate (Rate Gryoscope)

One of the most important measurements is the roll rate. The $x$ direction in our sensor frame lines up with the vertical axis of the rocket. So $x$ acceleration will be "Up" and gyro will be "roll":



In [2]:

    
fig = utils.Plot(r"Launch 11 Roll Rate", "Time [$s$]", "Rate [${}^0/s$]")
fig.plot(timestamp, gyr_x)
fig.ylim([-40, 425])
fig.show()

Oops! Unfortunately we saturated our roll sensor. But should be able to get some data out of the magnetometer.

Since the magetometer is looking mostly at the Earth's relativly fixed magnetic field, I would expect to see it change as the rocket spins. Here we want $y$ and $z$ axis:



In [3]:

    
fig = utils.Plot(r"Launch 11 Magnetometer", "Time [$s$]", "Magnetometer Values [$T$]")
fig.plot(timestamp, mag_y, label=r"$y$")
fig.plot(timestamp, mag_z, color=utils.blue, label=r"$z$")
fig.legend(loc=3)
fig.show()

If we treat this pair of numbers as a vector, then we should see the vector run around the origin as the rocket spins:



In [4]:

    
# Trick! Get the orgin to be zero by setting the mean to be zero
mag_y = mag_y-mag_y.mean()
mag_z = mag_z-mag_z.mean()

fig = utils.Plot(r"Launch 11 Magnetometer: $[y,z]$ vector", "Time [$s$]", "Magnetometer [$T$]")
fig.plot(mag_y, mag_z)
fig.plot(mag_y[:100], mag_z[:100], color='#00ff00') # mark the begining, look for a small green dot
fig.show()

Yep! That's some spinnin. Now we can compute roll angle:



In [5]:

    
# use numpy to build [y,z] vector from the two arrays
mag = np.dstack((mag_y, mag_z))[0]

# reference angle, should be a mean, but good enough.
begin = mag[0]

# storage
angles = []     # roll angle
rollup = 0      # every time we make one revolution, add 360 to this
last_angle = 0  # for checking if we've wrapped back to -pi

for v1 in mag:
    
    # compute angle between now and the reference angle
    s = np.cross(v1, begin)
    c = np.dot(v1, begin)
    angle = np.arctan2(s, c)
    
    # if we just wrapped around, add a revolution
    if last_angle > 2 and angle < 0:
        rollup += 360
        
    # store
    angles.append(np.degrees(angle)+rollup)
    last_angle = angle
    
# how many times did we spin?
print "The rocket made", rollup/360, "revolutions during ascent"

fig = utils.Plot(r"Launch 11 Magetometer Derived Roll Angle", "Time [$s$]", "Total Angle [${}^0$]")
fig.plot(timestamp, angles)
fig.ylim([-200,14000])
fig.show()









    



The rocket made 38 revolutions during ascent

We can differentiate this to get the aproximate roll rate. But! It will be very noisy



In [6]:

    
# differentiate
sample_rate = 819.2
angle_rate = np.diff(angles)
angle_rate = angle_rate*sample_rate  # correct for time units using the sample rate of the sensor

# low pass the hell out of the result
from scipy.signal import firwin, lfilter
freq_cuttoff = 0.5
filter_taps = 4096
window = firwin(numtaps=filter_taps, cutoff=freq_cuttoff, nyq=sample_rate/2)
angle_rate_padded = np.concatenate((angle_rate, np.zeros(filter_taps / 2)))
angle_rate_filtered = lfilter(window, 1.0, angle_rate_padded)[filter_taps / 2:]


fig = utils.Plot(r"Launch 11 Magetometer Derived Roll Rate", "Time [$s$]", "Angle [${}^0/s$]")
fig.plot(timestamp[1:], angle_rate, color='black', alpha=0.1, label="Roll Rate")
fig.plot(timestamp[1:], angle_rate_filtered, alpha=1, lw=1.6, label="Filtered Roll Rate")
fig.plot(timestamp, gyr_x, color=utils.blue, label="Roll Gryo")
fig.ylim([-80, 1000])
fig.legend()
fig.show()

Unfortunatly (again) after so much filtering, none of the subtle effects will show up in the filtered magnetometer data. But it's interesting to see the total roll rate, which orgininally I thought was lost to sensor saturation.