Canard Aerodynamics

We have 4 small "delta" wing canards. We want to model the forces on them, specifically lift.

Lift

The amount of lift on a wing can be very simply modeled by the following equation:

$$\begin{equation} L = C_{L\alpha}\frac{1}{2}\rho v^2 A \end{equation}$$

Where

  • $L$ is the lift
  • $C_{L\alpha}$ is the coefficient of lift (at a given angle of attack)
  • $\rho$ is the density of the atmosphere
  • $v$ is the velocity of the rocket
  • $A$ is the planform area of the fin

All of the complexity is hidden in the $C_{L\alpha}$ term! In fact this equation is a little disingenuous. What happened is $\frac{1}{2}\rho v^2 A$ gives the wrong answer for lift, so we conveniently pick a $C_{L}$ so that we magically get the right lift!

This is pretty useful when looking at wind tunnel data. You can trivially compute a lift coefficient from the measured air speed, air density, wing size, and lift. Save your measurements to a lookup-table and tada! You can get the lift anytime you want.

We don't have a windtunnel. So how are we going to get reasonable numbers for $C_{L}$? Much has been written on the subject. There are no analytical soulutions for general shapes. But we have a very specific shape in mind. For our kind of wing wind tunnel data has been taken and condensed into analytical aproximations over a wide range of speeds and angles of attack.

Subsonic

Lets start with the subsonic case (< Mach 0.8). Lift coefficients are very sensitive to angle of attack, so we'll see that as the main dependency.

$$\begin{equation} C_{L\alpha} = K_P\sin\alpha\cos^2\alpha + K_V\cos\alpha\sin^2\alpha \end{equation}$$

Where $K_P$ and $K_V$ are more magic constants. However these constants can are related to the fin geometry, and don't change during flight. Also note that $C_{L\alpha}$ will go to 0 at $\alpha$ = 0. Convinient!

Based on a lookup-table from published resutls, we're guessing the constants are:

  • $K_P$ = 2.45
  • $K_V$ = 3.2

Filling in those constants will let us now easily compute lift at low speeds!

Supersonic

The supersonic is a different regime. The aproximation now switches to being linear with angle of attack (really because of small-angle theorem). So let's define a new variable $C_{L_{base}}$ that is defined as:

$$\begin{equation} C_{L_{base}} \equiv \frac{ C_{L\alpha} }{\alpha} \end{equation}$$

So now all we need to determine is $C_{L_{base}}$!

...

Apparently it's about 3.2 at around Mach 1. We don't spend much time in the supersonic, and we don't get much above Mach 1, so for our purposes this is probably enough.

Lift Coefficient

We should now have everything we need, so lets look at $C_L$


In [1]:
from math import sin, cos, radians, exp

# Define PSAS wing:
k_p = 2.45
k_v = 3.2
Cl_base = 3.2
fin_area = 1.13e-3


def C_L(a, v):
    """Find C_L for a given speed and angle of attack
    
    :param float a: Angle of attack, alpha in degrees
    :param float v: velocity, v in m/s
    """
    
    # math is in radians
    a = radians(a)
    
    # Subsonic case
    def _subsonic():
        sina = sin(a)
        cosa = cos(a)
        cl = k_p*sina*cosa*cosa
        cl += k_v*cosa*sina*sina
        return cl
    
    # Supersonic case
    def _supersonic():
        cl = a*Cl_base
        return cl

    if v <= 265:
        return _subsonic()
    elif v < 330:
        # Intepolate between super and subsonic
        y0 = _subsonic()
        y1 = _supersonic()
        x0 = 265
        x1 = 330
        cl = y0 + (y1-y0)*(v - x0)/(x1-x0)
        return cl
    else:
        return _supersonic()

Now lets graph the results over a couple of cases:


In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams.update({'font.size': 16})

In [3]:
# range of angle of attack to chart
alpha = [a*0.5 for a in range(31)]
CL = [C_L(a, 130) for a in alpha]

fig = plt.figure(figsize=(15,9))
ax = fig.add_subplot(111)
plt.title("Angle of Attack effect on CL at Mach=0.4")
plt.ylabel(r"$C_{L\alpha}$")
plt.xlabel(r"$\alpha$ [${}^0$]")

plt.plot(alpha, CL)

ax.grid(color='grey', alpha=0.3, linestyle='dashed', linewidth=0.5)
plt.show()



In [4]:
# range of angle of attack to chart
alpha = [a*0.5 for a in range(31)]
CL = [C_L(a, 360) for a in alpha]

fig = plt.figure(figsize=(15,9))
ax = fig.add_subplot(111)
plt.title("Angle of Attack effect on CL at Mach=1.1")
plt.ylabel(r"$C_{L\alpha}$")
plt.xlabel(r"$\alpha$ [${}^0$]")

plt.plot(alpha, CL)

ax.grid(color='grey', alpha=0.3, linestyle='dashed', linewidth=0.5)
plt.show()



In [5]:
# Range of mach's
vel = [0.01*m for m in range(70,110,1)]
CL = [C_L(5, v*330) for v in vel]

fig = plt.figure(figsize=(15,9))
ax = fig.add_subplot(111)
plt.title(r"Mach effect on $C_L$ at $\alpha = 5^0$")
plt.ylabel(r"$C_{L\alpha}$, $\alpha = 5^0$")
plt.xlabel(r"Mach number")

plt.plot(vel, CL, 'r.-')

ax.grid(color='grey', alpha=0.3, linestyle='dashed', linewidth=0.5)
plt.show()


Lift!

Now we have our approximation for lift on one canard! We can lookup $C_L$ useing the above C_L() function.


In [6]:
def lift(a, v, alt):
    """Compute the lift of one fin at an angle of
    attack, velocity, and altitdue
    
    :param float a: Angle of attack in degrees
    :param float v: velocity in m/s
    :param float alt: altitude MSL in meters
    :returns float lift in Newtons:
    
    """
    
    # get density of atmosphere with quick exponential model
    rho = 1.2250 * exp((-9.80665 * 0.0289644 * alt)/(8.31432*288.15))

    l = 0.5*C_L(a, v)*rho*v*v*fin_area
    
    return l

Example Lift Curves


In [7]:
# Lift as a function of velocity at a=5deg
vel = [0.01*m for m in range(10,120,1)]
l = [lift(5, v*330, 3000) for v in vel]

fig = plt.figure(figsize=(15,9))
ax = fig.add_subplot(111)
plt.title(r"Lift vs Velocity at $\alpha = 5^0$")
plt.ylabel(r"Lift [N]")
plt.xlabel(r"Mach number")

plt.plot(vel, l, 'r.-')

ax.grid(color='grey', alpha=0.2, linestyle='dashed', linewidth=0.5)
plt.show()



In [8]:
vel = [0.01*m for m in range(10,120,1)]
l = [lift(1, v*330, 3000) for v in vel]

fig = plt.figure(figsize=(15,9))
ax = fig.add_subplot(111)
plt.title(r"Lift vs Velocity at $\alpha = 1^0$")
plt.ylabel(r"Lift [N]")
plt.xlabel(r"Mach number")

plt.plot(vel, l, 'r.-')

ax.grid(color='grey', alpha=0.2, linestyle='dashed', linewidth=0.5)
plt.show()