We have 4 small "delta" wing canards. We want to model the forces on them, specifically 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
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.
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:
Filling in those constants will let us now easily compute lift at low speeds!
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.
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()
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
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()