Please note that we are stealing this aero model wholesale from OpenRocket with very little understanding of what is going on. Please read chapter 3 of the OpenRocket technical documentation.
This is the Barrowman method, as extended by Galejs.
Angle of attack: displacement of rocket's z-axis from relative velocity vector
Normal Force: Corrective force normal to z-axis of rocket
In general, the normal and drag forces are not orthogonal. They are only when angle of attack is 0. We define two other forces to produce distinct pairs of orthogonal forces.
Side force: Orthogonal to drag, normal to direction of motion
Roll: around z-axis of rocket
A consequence of Barrowman's thesis is that cylindrical body components do not effect normal forces or the center of pressure location.Also, for the cylinder, the cross-sectional area is not changing w.r.t. axial position, so there is no contribution to pitching moment either. However, an extension to the method shows that cylinders contribute to body lift forces when angle of attack is not small (i.e., ignored).
We have at least two fins and they aren't canted, so we have no yawing or rolling moments from the fins.
We ignore fins interference with body but calculate body's interference with fins.
We ignore pitch damping torque since it is only relevant at apogee.
Our fins are not canted, so there is no roll forcing torque.
We have done most of our calculations for Normal Forces, which are paired with Axial Drag. It will be useful for simulation purposes to break Normal Forces into components and calculate Drag forces as well which are distinct, and describe deceleration.
Barrowman gives methods for both fully turbulent and partially-laminar layers, but OpenRocket assumes boundary layers are fully turbulent (due to < 5% difference in apogee).
In [1]:
import scipy.interpolate
import numpy as np
class AeroModel():
def __init__(self, diameter, body_len, nose_len, fin):
# kinematic viscosity of air
self.mu = 1.5 * 10**-5 # m^2/s
#self.R_crit = 5*10**5 # critical reynolds number at which flow becomes turbulent
#self.sincAOA = 1.0 # to avoid avoiding dividing by zero, if we need it
self.painted_surface = 2 * 10**-6 # m
self.reference_len = diameter
self.body_len = body_len
self.nose_len = nose_len
self.A_ref = np.pi* self.reference_len**2 / 4
self.R_crit_skin = 51 * (self.painted_surface/self.body_len)**-1.039
self.height = fin.height # fin span
self.root = fin.root
self.fin_lead = body_len - self.root
self.tip = fin.tip
self.thickness = fin.thickness
self.sweep_angle = fin.sweep
self.fin_area = fin.volume/self.thickness
self.MAC_length = 2/3 * (self.root**2 + self.tip**2 + self.root*self.tip)/(self.root + self.tip)
self.MACspan = self.MAC_span()
self.MAC_lead = self.effective_MAC_LE()
self.finAR = self.fin_aspect()
# Commercial Airplane Design Principles, eq 5.11 for midchord sweep
self.midchord_sweep = np.arctan(self.sweep_angle - (2/self.finAR) * (1-self.tip/self.root)/(1+self.tip/self.root))
self.fin_poly = self.coefficients(self.finAR)
self.fin_pressure_mult = np.cos(self.sweep_angle)**2 * 4 * self.thickness * self.height / self.A_ref
self.bodytail = self.bodytail_interference()
self.finshape = self.roll_damp_fin_shape()
self.fin_interps = self.Kvalues()
self.CDfric_coef = self.frictcoef()
def fin_aspect(self):
return 2*self.height**2/self.fin_area # multiplying by 2 because openrocket does, technically that may be wrong
#fin-body interference coefficient
def bodytail_interference(self):
return 1 + (self.reference_len/2)/(self.height + self.reference_len/2)
# subsonic fin shape term, eq 3.70
def roll_damp_fin_shape(self):
radius = self.reference_len/2
term1 = (self.root + self.tip) * self.height *radius**2 / 2
term2 = (self.root + 2*self.tip) * self.height**2 * radius/3
term3 = (self.root + 3*self.tip) * self.height**3 / 12
return term1+term2+term3
# skin friction drag coefficient, assuming nose cone is cylinder for now
# func of body fineness, fin thickness, mean aero chord, wet area of body, wet area of fins (both sides incl)
def frictcoef(self):
return ((1 + 1/(2*(self.body_len/self.reference_len)))* (np.pi * self.reference_len * self.body_len) +
(1 + 2*self.thickness/self.MAC_length)* ((2*self.fin_area+self.tip*self.thickness)*4)) / self.A_ref
# reciprocal of Prandtl factor P "which corrects subsonic force coefficients for compressible flow"
def Beta(self, mach):
return np.sqrt(abs(mach**2 - 1))
# "v0 is the free-stream velocity of the rocket, ν is the kinematic viscosity of air"
def Reynold(self, v0):
return self.body_len*v0/self.mu
#fin interpolator, yikes copypasta
def coefficients(self, ar):
denom = (1- 3.4641*ar)**2
poly=[]
poly.append((9.16049 * (-0.588838 + ar) * (-0.20624 + ar)) / denom)
poly.append((-31.6049 * (-0.705375 + ar) * (-0.198476 + ar)) / denom)
poly.append((55.3086 * (-0.711482 + ar) * (-0.196772 + ar)) / denom)
poly.append((-39.5062 * (-0.72074 + ar) * (-0.194245 + ar)) / denom)
poly.append((12.8395 * (-0.725688 + ar) * (-0.19292 + ar)) / denom)
poly.append((-1.58025 * (-0.728769 + ar) * (-0.192105 + ar)) / denom)
return poly
# coefficients for single fin supersonic normal force coefficient derivative
# we will follow their approach and precompute some coefficients
# and we will return linear interpolators
def Kvalues(self):
gamma = 1.4 # wtf is this
Ma0 = 1.5
def K1(beta):
return 2.0/beta
def K2(beta, mach):
return ((gamma + 1) * mach**4 - 4 * beta**2) / (4 * beta**4)
def K3(beta, mach):
return ((gamma + 1) * mach**8 + (
2 * gamma**2 - 7 * gamma - 5) * mach**6
+ 10 * (gamma + 1) * mach**4 + 8) / (6 * beta**7)
Ma, k1, k2, k3 = [Ma0], [], [], []
while Ma[-1] < 5.0:
beta = self.Beta(Ma[-1])
k1.append(K1(beta))
k2.append(K2(beta, Ma[-1]))
k3.append(K3(beta, Ma[-1]))
Ma.append(Ma[-1] + 0.1)
trash = Ma.pop()
k1_intrp = lambda M: np.interp(M, Ma, k1)
k2_intrp = lambda M: np.interp(M, Ma, k2)
k3_intrp = lambda M: np.interp(M, Ma, k3)
# sneaking a second interpolator creation in here
self.CD_interp_angle = 17 * np.pi / 180
start = (0,1)
end = (self.CD_interp_angle, 1.3)
self.CD_interp_below = scipy.interpolate.CubicHermiteSpline([start[0], end[0]],
[start[1], end[1]],
[0, 0],
extrapolate=False)
start = (self.CD_interp_angle, 1.3)
end = (np.pi/2, 0)
self.CD_interp_above = scipy.interpolate.CubicHermiteSpline([start[0], end[0]],
[start[1], end[1]],
[0, 0],
extrapolate=False)
# sneaking some values for a third interpolator here
self.subsonic = 0.9
self.supersonic = 1.5
self.supersonic_b = (self.supersonic**2 - 1)**1.5
self.sq = np.sqrt(1 +
(1 - self.subsonic**2) * (self.height**2 /
(self.fin_area * np.cos(self.midchord_sweep)))**2)
self.subV = 2 * np.pi * self.height**2 / self.A_ref / (1 + self.sq)
self.subD_coeff = 2 * np.pi * self.height**6 / ((self.fin_area * np.cos(self.midchord_sweep))**2
* self.A_ref * self.sq * (1 + self.sq)**2)
self.superD = - self.fin_area/self.A_ref * 2 * self.supersonic / self.supersonic_b
self.superV = None
# sneaking in a 4th interpolator
self.CN_a1_sub = lambda b: 2*np.pi*self.height**2/self.A_ref/(1+np.sqrt(
1 + (self.height**2*b/(self.fin_area*np.cos(self.midchord_sweep)))**2))
# oh, a 5th and 6th interpolator. for von Karman/LD-Haack nosecone pressure drag as function of mach
# https://github.com/openrocket/openrocket/blob/unstable/core/src/net/sf/openrocket/aerodynamics/barrowman/SymmetricComponentCalc.java
self.nose_extrap = np.log(self.nose_len / self.A_ref + 1) / np.log(4)
self.vonKarmanInterpolator = scipy.interpolate.interp1d(
[0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4, 1.6, 2.0, 3.0],
[0, 0.010, 0.027, 0.055, 0.070, 0.081, 0.095, 0.097, 0.091, 0.083],
kind='linear', copy=False, fill_value="extrapolate")
# this block unnessary, just look at the first value
# only leaving it for now in case we wanna different nosecone
minValue = self.vonKarmanInterpolator.__call__(0.9)
if minValue > 0.001:
minDeriv = (self.vonKarmanInterpolator.__call__(0.9 + 0.01) - minValue) / 0.01
a = minValue - 0.8 * np.sin(0.00)**2
b = minDeriv / a
self.nose_cone_pressure_subsonic = lambda M: a * M**b + 0.8 * np.sin(0.00)**2
else:
self.nose_cone_pressure_subsonic = lambda M: 0
return (k1_intrp, k1_intrp, k1_intrp)
# assuming CP constant, fin for mach < 0.5
# for trapezoid fins with tip parallel to root
def MAC_span(self):
return (self.height/3) * (self.root + 2*self.tip)/(self.root + self.tip)
def effective_MAC_LE(self):
def Z_t():
return np.tan(self.sweep_angle)*self.height
return (Z_t()/3)*(self.root + 2*self.tip)/(self.root + self.tip)
# correction for fin CP changing with speed
# for mach above 2, and below is interpolated
# by a 5th order overfit polynomial lol. RHS = right hand side of equation
def fin_CoP(self, mach):
if mach <=0.7: #subsonic, 1/4 of MAC
RHS = 0.25
# return self.fin_lead + self.MAC_lead + self.MAC_length/4
elif mach >= 2: #supersonic, 1/2 of mac
beta = self.Beta(mach)
RHS = (self.finAR*beta-0.67)/(2*self.finAR*beta-1)
else: #transonic, oh god
x=1.0
RHS=0
for coeff in self.fin_poly:
RHS += coeff*x
x *= mach
return self.fin_lead + self.MAC_lead + self.MAC_length*RHS
#Tactical Missile Design by Fleeman, ignore variance with mach
# this is from slender body and crossflow theories # deprecated?
def bodynose_CoP(self, alpha):
sin2a = np.sin(alpha)**2
return 0.63*self.nose_len*(1-sin2a) + 0.5*self.body_len*sin2a
# robert galejs' correction term for cylindrical body section lift
# reference area is cross-section, planform is diam*length
# actually this comes from pg 3-11 of Fluid Dynamic Drag, Hoerner
def C_N_alpha_lift(self, alpha):
K = 1.1 # fudge factor
return K*(self.body_len-self.nose_len)*self.reference_len/self.A_ref * np.sin(alpha) * np.sinc(alpha/np.pi)
# nosecone coefficients
def C_N_alpha_nose(self, alpha):
return 2*np.sinc(alpha/np.pi)
def normal_force_coefficient(self, alpha, mach):
# N=4 fins equally spaced, N/2
def C_N_alpha_fins(alpha, mach):
# single fin subsonic normal force coefficient derivative
def C_N_alpha_1_sub(beta):
return self.CN_a1_sub(beta)
# single fin supersonic normal force coefficient derivative
# Barrowman used a third-order (!) expansion of Busemann theory
# openrocket simplified his method by assuming fins are flat plates
# and ignoring correction of fin-tip mach cones
def C_N_alpha_1_super(alpha, mach):
return self.fin_area/self.A_ref * (self.fin_interps[0](mach)
+ self.fin_interps[1](mach) * alpha
+ self.fin_interps[2](mach) * alpha**2)
# black magic, don't ask
def C_N_alpha_1_trans(alpha, mach):
if not self.superV:
self.superV = lambda a: C_N_alpha_1_super(a, self.supersonic)
trans_interp = scipy.interpolate.CubicHermiteSpline([self.subsonic, self.supersonic],
[self.subV, self.superV(alpha)],
[self.subD_coeff*mach, self.superD])
return trans_interp.__call__(mach)
if mach <= 0.9:
CNa1 = C_N_alpha_1_sub(self.Beta(mach))
elif mach > 1.5:
CNa1 = C_N_alpha_1_super(alpha, mach)
else: #cry
CNa1 = C_N_alpha_1_trans(alpha, mach)
return CNa1 * 2
#fin-body interference
def C_N_alpha_tail(CNa_fins):
return CNa_fins * self.bodytail
return C_N_alpha_tail(C_N_alpha_fins(alpha, mach))
# deprecated
def C_m_alpha_nose(self, alpha): # for now, not parametrized
return 2/(self.reference_len*self.A_ref) * (0.5*self.A_ref -
np.pi/3 * 0.5 * (self.reference_len/2)**2) * np.sinc(alpha/np.pi)
def roll_damping_coefficient(self, omega, v0, Mach):
# subsonic roll damping from 4 fins, eq 3.69
# v0 is free-stream velocity, omega is roll rate
def C_ld_subsonic(Mach, v0, omega):
CNa0 = 2*np.pi / self.Beta(Mach)
v0 = v0 if v0 != 0 else 1
return 4*CNa0*omega*self.finshape/(self.A_ref*self.reference_len*v0)
def C_ld_supersonic(mach, omega, v0):
# local pressure coefficient of a strip of fin
def C_P_i(eta, K):
return (K[0] * eta
+ K[1] * eta**2
+ K[2] * eta**3)
def eta_i(distance, coeff):# eq 3.62, small angle approx
return distance*coeff # reasonable up to 20 deg
def width(y): #trapezoid, y is height
return y*self.tip/self.height + (self.height-y)*self.root/self.height
vel_coeff = omega/v0 if v0 != 0 else omega
eval_K = [k(mach) for k in self.fin_interps]
mult = 4/(self.A_ref*self.reference_len)
num_terms = 5
del_y = self.height/num_terms
terms = sum([C_P_i(eta_i(
i*del_y, vel_coeff),
eval_K) *i*del_y * width(i*del_y) * del_y
for i in range(num_terms+1)])
return mult*terms
if omega == 0: return 0
if Mach <= 1:
return C_ld_subsonic(Mach, v0, omega)
else:
return C_ld_supersonic(Mach, omega, v0)
def axial_drag_coefficient(self, R, M, alpha):
# drag coefficient at zero angle of attack
def C_D_0(R, M):
# skin friction coefficient function of reynolds number and mach
# turbulent flow, in BarrowmanCalculator.java there is a continuous interpolation of correction around Mach 1
# but that probably isn't needed here. could precalc a little here?
def eh_C_f(R, M):
#if M <= 1:# compressibility correction
# correction = 1 - 0.1*M**2
#elif M > 1:
# correction = (1 + 0.15 * M**2)**-0.58
if M < 0.9:# compressibility correction
correction = 1 - 0.1*M**2
elif M > 1.1:
correction = 1/(1 + 0.18 * M**2)
else:
c1 = 1 - 0.1*M**2
c2 = 1/(1 + 0.18 * M**2)
correction = c2 * (M - 0.9) / 0.2 + c1 * (1.1 - M) / 0.2
if R < 10**4: # velocity < 1 m/s
return 1.48 * 10**-2
elif R < self.R_crit_skin: #assuming smooth surface
return correction * (1.5* np.log(R) - 5.6)**-2
else: # simplified form of eqs (3.80) and (3.79)
if M > 1: correction = max(1/(1+ 0.18 * M**2), correction)
return correction * 0.032*(self.painted_surface / self.body_len)**0.2
#return correction * 0.06821 * self.R_crit_skin**-.1925
def C_f(R, M):
correction1, correction2 = 1, 1
if M < 1.1:# compressibility correction
correction1 = 1 - 0.1*M**2
if M > 0.9:
correction2 = (1 + 0.15 * M**2)**-0.58
if M < 0.9:
correction = correction1
elif M < 1.1:
correction = correction2*(M-.9)/.2 + correction1*(1.1-M)/.2
else:
correction = correction2
if M < 0.9:# compressibility correction
roughcorrection = 1 - 0.1*M**2
elif M > 1.1:
roughcorrection = 1/(1 + 0.18 * M**2)
else:
c1 = 1 - 0.1*M**2
c2 = 1/(1 + 0.18 * M**2)
roughcorrection = c2 * (M - 0.9) / 0.2 + c1 * (1.1 - M) / 0.2
if R < 10**4: # velocity < 1 m/s
Cf = 1.48 * 10**-2
else: #assuming smooth surface
Cf = correction * (1.5* np.log(R) - 5.6)**-2
roughlimited = roughcorrection*0.032*(self.painted_surface / self.body_len)**0.2
Cf = max(Cf, roughlimited)
return Cf
#else: # simplified form of eqs (3.80) and (3.79)
# if M > 1: correction = max(1/(1+ 0.18 * M**2), correction)
# return correction *
#return correction * 0.06821 * self.R_crit_skin**-.1925
def new_C_f(R, M):
if M < 0.9:# compressibility correction
correction = 1 - 0.1*M**2
elif M > 1.1:
correction = 1/(1 + 0.18 * M**2)
else:
c1 = 1 - 0.1*M**2
c2 = 1/(1 + 0.18 * M**2)
correction = c2 * (M - 0.9) / 0.2 + c1 * (1.1 - M) / 0.2
if R < 10**4: # velocity < 1 m/s
return 1.33 * 10**-2
elif R < 5.39 * 10**5: #assuming smooth surface
return correction * 1.328 / np.sqrt(R)
else: # simplified form of eqs (3.80) and (3.79)
if M > 1: correction = max(1/(1+ 0.18 * M**2), correction)
return correction * ((1.5* np.log(R) - 5.6)**-2 - 1700/R)
#return correction * 0.06821 * self.R_crit_skin**-.1925
# skin friction drag coefficient
# func of body fineness, fin thickness, mean aero chord, wet area of body, wet area of fins (both sides incl)
def C_D_frict(Cfc):
return Cfc * self.CDfric_coef
# nosecone pressure drag coefficient
# assumes nosecone smoothly transitions
def C_D_nose_press(M):
def stag(M): # eq B.1 and B.2, blunt cylinder stagnation pressure
if M <= 1:
return 0.85 * (1 + M**2 / 4 + M**4 / 40)
else:
return 0.85 * (1.84 - 0.76 / M**2 + 0.166 / M**4 + 0.035 / M**6)
if M < 0.9:
# temporary hack
#C0 = stag(.9)
#minValue = C0 * (self.vonKarmanInterpolator.__call__(0.9) / C0)**self.nose_extrap
#minDeriv = (C0 * (self.vonKarmanInterpolator.__call__(0.9 + .01) / C0)**self.nose_extrap - minValue) / 0.01
#b = minDeriv / minValue
#temporary = lambda M: minValue * M**b
#return temporary(M)
return self.nose_cone_pressure_subsonic(M) # this should be right
else:
C0 = stag(M)
# ref eq B.9
return C0 * (self.vonKarmanInterpolator.__call__(M) / C0)**self.nose_extrap
#return .1005025 + (.002503454 - .1005025)/(1+ (M/1.065884)**6.060682)
# fin pressure drag coefficient, eq 3.89, 3.91
# assumes leading edge rounded and trailing edge tapered
def C_D_fin_press(M):
if M < 0.9:
LE = (1 - M**2)**-0.417 - 1
elif M < 1:
LE = 1 - 1.785*(M - 0.9)
else:
LE = 1.214 - 0.502 / M**2 + 0.1095 / M**4
return LE * self.fin_pressure_mult # later check ref areas again
# base drag coefficient
# eq 3.94, at some point need to think about whether this only applies to coasting phase
def C_D_base(M):
if M < 1:
return 0.12 + 0.13*M**2
else:
return 0.25/M
# note we are ignoring parasitic drag from launch lugs/pipes for now
Cfc = C_f(R, M)
CD_f = C_D_frict(Cfc)
CD_p1 = C_D_nose_press(M)
CD_p2 = C_D_fin_press(M)
CD_b = C_D_base(M)
return CD_f + CD_p1 + CD_p2 + CD_b
# drag coefficient adjusted for angle of attack
def C_D_axial(alpha, CD0):
if alpha < 0.01: return CD0
alpha_ = np.pi - alpha if alpha > np.pi/2 else alpha
if alpha_ <= self.CD_interp_angle:
mult = self.CD_interp_below.__call__(alpha_)
else:
mult = self.CD_interp_above.__call__(alpha_)
if alpha < np.pi/2:
return mult * CD0
else:
return -mult * CD0
return C_D_axial(alpha, C_D_0(R, M))
def physics(self, alpha, Mach, v0, omega, mu):
self.mu = mu
CoPs = [self.bodynose_CoP(alpha), self.fin_CoP(Mach)]
coeffs = [self.C_N_alpha_lift(alpha)+ self.C_N_alpha_nose(alpha), self.normal_force_coefficient(alpha, Mach)]
CoP = sum([CoPs[i]* CN_alpha for i, CN_alpha in enumerate(coeffs)]) / sum(coeffs)
# forces
CN = alpha *sum(coeffs)
CDax = self.axial_drag_coefficient(self.Reynold(v0), Mach, alpha)
# moments
Cld_times_d = self.roll_damping_coefficient(omega, v0, Mach) * self.reference_len
return CoP, CN, CDax, Cld_times_d # remember that CoP is referenced from nosecone, not from engine
In [ ]: