Okay, for the most part, this isn't analyzing drag for the recovery system yet, lol. Most of it is looking at the OpenRocket outputs and comparing them to DATCOM. In fact, this will probably get blobbed together as some omni-analysis for LV3 at some point.
In [1]:
gamma=1.4 # C_p/C_v for air, taken to be constant
In [2]:
getwd()
if (!exists('par.default'))
par.default <- par(no.readonly=T)
openRocketCSV <- "simData/LV3_L13_ideal_2017-5-30-1804.csv"
dat <- read.csv(openRocketCSV, skip=6)
datNames <- names(dat)
dat <- read.table(openRocketCSV, header= F, sep=',')
names(dat) <- datNames
head(dat)
# Print out the event lines from the CSV:
# system('grep \'occur\' simData/LV3_L13_ideal_2017-5-30-1804.csv', intern = TRUE)
events <- as.numeric(system('grep \'occur\' simData/LV3_L13_ideal_2017-5-30-1804.csv | grep -o \'[0-9]\\+\\.*[0-9]*\'', intern = TRUE))
events.names <- system('grep \'occur\' simData/LV3_L13_ideal_2017-5-30-1804.csv | grep -o \'[[:upper:]_]\\{2,\\}\'', intern = TRUE)
names(events) <- events.names
events
In [3]:
# Print out the event lines from the CSV:
# system('grep \'occur\' simData/LV3_L13_ideal_2017-5-30-1804.csv', intern = TRUE)
events <- as.numeric(system('grep \'occur\' simData/LV3_L13_ideal_2017-5-30-1804.csv | grep -o \'[0-9]\\+\\.*[0-9]*\'', intern = TRUE))
events.names <- system('grep \'occur\' simData/LV3_L13_ideal_2017-5-30-1804.csv | grep -o \'[[:upper:]_]\\{2,\\}\'', intern = TRUE)
names(events) <- events.names
# events
In [4]:
dat$colors <- sapply(dat$X..Time..s., function(t) sum(t>=events))
# layout(mat = matrix(c(1,2), nrow = 1))
plot(dat$X..Time..s., dat$Altitude..km., col=dat$colors, pch=20)
abline(v=events, col=1:length(events))
# text(
# labels=names(events),
# x=events,
# y=seq(from=0, to=max(dat$Alt), length.out=length(events)),
# col=1:length(events),
# pos=4, srt=0, offset=-0.0
# )
legend(
'right',
legend=paste(names(events), 'at t=', events),
col=1:length(events),
pch=20,
title='Colors highlight times following their event.'
)
In [5]:
ind <- which(
dat$X..Time..s. > events['LAUNCHROD']
& dat$X..Time..s. < events['APOGEE'] +2
)
launchlim <- range(dat$X..Time..s.[ind])
In [6]:
layout(matrix(c(1,2),nrow= 2))
plot(dat$X..Time..s, dat$Drag.coefficient, col=dat$col, xlim=launchlim)
abline(v=events, col=1:length(events))
plot(dat$X..Time..s., dat$Total.velocity..m.s., col=dat$col, xlim=launchlim)
abline(v=events, col=1:length(events))
Scheduled test points are shown as numbers over the plots.
In [7]:
rodAlt <- dat$Altitude..km.[which(dat$X..Time..s.==events['LAUNCHROD'])]
ejectionAlt <- dat$Altitude..km.[which(dat$X..Time..s.==events['EJECTION_CHARGE'])]
alt.schedule.input <- c(0,rodAlt,0.02,0.1,0.4,0.8,1.4,ejectionAlt,3,5,7,max(dat$Alt))
# find the first index where we reach that altitude:
ind.schedule <- sapply(alt.schedule.input, function(x) which.max(dat$Alt >= x))
In [8]:
# schedule <- data.frame(matrix(ncol=0, nrow=length(alt.schedule)))
# for (name in names(dat))
# {
# tempsch <- data.frame(approx(dat$Alt, dat[[name]], xout=alt.schedule, rule=2)$y)
# names(tempsch) <- name
# schedule <- cbind(schedule, tempsch)
# }
schedule <- dat[ind.schedule,]
schedule$test.number <- 1:length(ind.schedule)
schedule
write.csv(schedule, file='testPoints_OpenRocket.csv')
In [9]:
# rodAlt <- dat$Altitude..km.[which(dat$X..Time..s.==events['LAUNCHROD'])]
# ejectionAlt <- dat$Altitude..km.[which(dat$X..Time..s.==events['EJECTION_CHARGE'])]
# schedule.alt.mach <- approx( # sample some altitude/mach number pairs with linear interpolation
# dat$Alt[ind], dat$Mach[ind],
# xout=c(0,rodAlt,0.02,0.1,0.4,0.8,1.4,ejectionAlt,3,5,7,max(dat$Alt)), # hand-picked
# rule=2
# )
# alt.schedule <- schedule.alt.mach$x
# mach.schedule <- schedule.alt.mach$y
# # stability margin caliber schedule:
# smc.schedule <- approx(dat$Alt[ind], dat$Stability.margin.calibers[ind], xout=alt.schedule, rule=2)$y
# aoa.schedule <- approx(dat$Alt[ind], dat$Angle[ind], xout=alt.schedule, rule=2)$y
In [10]:
layout(matrix(c(1,2,3,4), nrow=4))
par.old <- par(mar=c(0,4,4,2)) # margins for top plot (bottom, left, top, right)
plot(dat$X..Time..s., dat$Altitude..km., xlim=launchlim, col=dat$col)
text(schedule$X..Time..s., schedule$Alt, schedule$test)
par(mar=c(0,4,0,2)) # margins for middle plot
plot(dat$X..Time..s., dat$Mach.number, xlim=launchlim, col=dat$col)
text(schedule$X..Time..s., schedule$Mach, schedule$test)
plot(dat$X..Time..s., dat$Stability.margin.calibers, xlim=launchlim, col=dat$col)
text(schedule$X..Time..s., schedule$Stability.margin, schedule$test)
par(mar=c(5,4,0,2)) # margins for bottom plot
plot(
dat$X..Time..s.[ind], dat$Angle.of.attack[ind],
xlim=launchlim,
ylim=c(0,4),
col=dat$col
)
text(schedule$X..Time..s., schedule$Angle, schedule$test)
abline(h=c(0,1,2.5,4.5), col='lightgray', lty=2)
cat('max angle of attack during normal flight (launch rod to deployment):', max(na.omit(dat$Angle.of.attack[ind])))
In [11]:
layout(matrix(c(1,2,3), nrow=3))
par(mar=c(0,4,4,2))
plot(dat$Altitude..km.[ind], dat$Mach.number[ind], col=dat$col[ind])
text(as.character(1:length(schedule$Alt)), x=schedule$Alt, y=schedule$Mach)
par(mar=c(0,4,0,2)) # margins for middle plot
plot(
dat$Altitude..km.[ind], dat$Stability.margin[ind],
col=dat$col[ind], ylim=c(0, max(na.omit(dat$Stab)))
)
text(as.character(1:length(schedule$Alt)), x=schedule$Alt, y=schedule$Stability.margin)
par(mar=c(5,4,0,2))
plot(dat$Alt[ind], dat$Angle[ind], col=dat$col[ind])
text(as.character(1:length(schedule$Alt)), x=schedule$Alt, y=schedule$Angle)
par(par.old)
layout(1)
In [12]:
# DATCOM-friendly output:
alt.schedule.inch <- schedule$Alt*1e6/25.4 #convert from km into in
cat('##### INSERT VALUES THESE INTO THE DATCOM INPUT FILE: #####\n')
cat('NMACH=')
cat(sprintf('%1.1f', length(schedule$Mach))); cat(',')
cat('\nMACH=')
cat(toupper(format(schedule$Mach, scientific = T, digits=2)),sep=','); cat(',')
cat('\nNALT=')
cat(sprintf('%1.1f', length(schedule$Alt))); cat(',')
cat('\nALT=')
cat(toupper(format(alt.schedule.inch, scientific = T, digits=2)),sep=','); cat(',')
In [13]:
plot(dat$Alt, dat$Drag.coefficient, col=dat$col)
# drag.schedule <- approx(dat$Alt, dat$Drag.coeff, xout=alt.schedule, rule=2)$y
grid()
text(as.character(1:length(schedule$Alt)), x=schedule$Alt, y=schedule$Drag.coeff)
In [14]:
galdat.drag <- read.csv('plots/Gallagher1947_40trapezoid_bevel_CDpoints.txt', comment.char = '#')
galdat.lift <- read.csv('plots/Gallagher1947_40trapezoid_bevel_CLpoints.txt', comment.char = '#')
In [15]:
names(galdat.drag) <- c('alpha', 'CD')
names(galdat.lift) <- c('alpha', 'CL')
galdat.drag
galdat.lift
In [16]:
m.drag <- lm('CD ~ I(alpha^3) + I(alpha^2) + alpha', data = galdat.drag)
m.lift <- lm('CL ~ I(alpha^3) + alpha + 0', data = galdat.lift)
# m.lift <- lm('CL ~ alpha', data = galdat.lift)
summary(m.drag)
summary(m.lift)
confint(m.drag)
confint(m.lift)
galdrag <- function(alpha) predict(m.drag, data.frame(alpha=alpha))
gallift <- function(alpha) predict(m.lift, data.frame(alpha=alpha))
plot(
galdat.drag$alpha, galdat.drag$CD,
xlim=c(0,50), ylim=c(0,1.3),
main='40 degree beveled trapezoidal wing, M=1.55
Gallagher and Mueller, 1947
(cubic fits)',
xlab='alpha (degrees)', ylab=NA
)
lines(0:50, galdrag(0:50))
points(galdat.lift$alpha, galdat.lift$CL, col='blue')
lines(0:50, gallift(0:50), col='blue')
legend(
'bottomright',
legend=c('drag', 'lift'),
col=c('black', 'blue'),
lty=1
)
In [17]:
semispan.exposed= 6.42*25.4/1e3 # meters
chord.root= 18*25.4/1e3 # meters
chord.tip= 5*25.4/1e3 # meters
fin.area <- semispan.exposed*(chord.root + chord.tip)/2
fin.area
$D = \frac{1}{2} \rho U^2 A C_D$
$L = \frac{1}{2} \rho U^2 A C_L$
In [18]:
ind.super <- dat$Mach >= 1.2
mask.super <- sapply(ind.super, function(x) {if (x) x else NA})
dat$fin.CD <- galdrag(dat$Angle)*mask.super
dat$fin.CL <- galdrag(dat$Angle)*mask.super
dat$density <- gamma*dat$Air.pressure..Pa./dat$Speed.of.sound..m.s.^2
dat$fin.Drag <- 1/2*dat$density*dat$Total.velocity..m.s.^2*fin.area*dat$fin.CD
dat$fin.Lift <- 1/2*dat$density*dat$Total.velocity..m.s.^2*fin.area*dat$fin.CL
In [19]:
dat$fin.normalLoad <- dat$fin.Lift*cos(dat$Angle*2*pi/360) + dat$fin.Drag*sin(dat$Angle*2*pi/360)
dat$fin.tangentialLoad <- dat$fin.Lift*sin(dat$Angle*2*pi/360) - dat$fin.Drag*cos(dat$Angle*2*pi/360)
# plot(dat$fin.norm, type='l')
# plot(dat$fin.tang, type='l', col='blue')
In [20]:
fin.normalLoad.max= 450*4.448 # N (converted from pounds)
FS.fin.normal <- fin.normalLoad.max/max(dat$fin.norm, na.rm = T )
Neglecting tangential loading on the fins, the factor of safety for normal (out-of-plane) loading is:
In [21]:
FS.fin.normal