Notice. New forum software under development. It's going to miss a few functions and look a bit ugly for a while, but I'm working on it full time now as the old forum was too unstable. Couple days, all good. If you notice any issues, please contact me.
|
Forum Index : Microcontroller and PC projects : MMX - model rocket performance
Author | Message | ||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
This post is a short computer program for the MicroMite eXtreme that estimates the flight performance of a single stage model rocket. Here's a typical user interaction with the software and the calculated results from the program. program rocket1 --------------- please input the launch site altitude (meters) (sites above sea level are positive, sites below sea level are negative) ? 0 please input the launch site temperature (degrees fahrenheit) ? 59 please input the thrust duration (seconds) ? 1.2 please input the total impulse (newton-seconds) ? 5 please input the initial mass (grams) ? 40 please input the propellant mass (grams) ? 8.33 please input the frontal diameter (millimeters) ? 18 please input the drag coefficient ? .321 program rocket1 --------------- burnout altitude 74.0670 meters burnout velocity 119.3594 meters per second burnout mass 31.6700 grams coast time 7.9316 seconds total flight time 9.1316 seconds maximum altitude 451.3929 meters average thrust 4.1667 newtons Here's the MMBASIC source code. ' program rocket1.bas March 19, 2017 ' MicroMite eXtreme version ' determines flight performance ' of single stage model rockets ' altitude at burnout, meters ' velocity at burnout, meters per second ' coast time and total flight time, seconds ' maximum altitude, meters ''''''''''''''''''''''''''' option default float dim altsite, tempsite, tduration, impulse, massi, mprop, fd, cd ' surface gravity (meters/second^2) const gravity = 9.80665 ' sea level atmospheric density (kilograms/cubic meter) const rhosl = 1.22557 print " " print "program rocket1" print "---------------" do print " " print "please input the launch site altitude (meters)" print "(sites above sea level are positive, sites below sea level are negative)" input altsite loop until (abs(altsite) >= 0.0) print " " print "please input the launch site temperature (degrees fahrenheit)" input tempsite do print " " print "please input the thrust duration (seconds)" input tduration loop until (tduration > 0.0) do print " " print "please input the total impulse (newton-seconds)" input impulse loop until (impulse > 0.0) do print " " print "please input the initial mass (grams)" input massi loop until (massi > 0.0) do print " " print "please input the propellant mass (grams)" input mprop loop until (mprop > 0.0) do print " " print "please input the frontal diameter (millimeters)" input fd loop until (fd > 0.0) do print " " print "please input the drag coefficient" input cd loop until (cd > 0.0) ' convert mass to kilograms and diameter to square meters massi = 0.001 * massi mprop = 0.001 * mprop fd = pi * fd * fd / 4000000.0 ' compensate for launch site altitude and temperature rho = rhosl * density(altsite) / (1.0 + (tempsite - 59.0) / 518.67) ' determine altitude performance thrust = impulse / tduration mass = (massi - 0.5 * mprop) k2 = 0.5 * rho * fd * cd weight = mass * gravity b = tduration * sqr(k2 * (thrust - weight)) / mass c = exp(b) d = exp(-b) e = 0.5 * (c + d) f = (c - d) / (c + d) ' calculate burnout conditions altbo = (mass / k2) * log(e) velbo = f * sqr((thrust - weight) / k2) massbo = massi - mprop weightbo = massbo * gravity ' compute coast conditions tcoast = sqr(massbo / (k2 * gravity)) * atn(velbo * sqr(k2 / weightbo)) altcoast = (massbo / (2.0 * k2)) * log(k2 * velbo * velbo / weightbo + 1.0) ' calculate total flight time and maximum altitude tflight = tduration + tcoast altmax = altbo + altcoast ' print results print " " print "program rocket1" print "---------------" print "burnout altitude ", str$(altbo, 0, 4), " meters" print " " print "burnout velocity ", STR$(velbo, 0, 4), " meters per second" print " " print "burnout mass ", STR$(1000.0 * massbo, 0, 4), " grams" print " " print "coast time ", str$(tcoast, 0, 4), " seconds" print " " print "total flight time ", str$(tflight, 0, 4), " seconds" print " " print "maximum altitude ", str$(altmax, 0, 4), " meters" print " " print "average thrust ", str$(thrust, 0, 4), " newtons" print " " end '''''''''''''''''''''''''' '''''''''''''''''''''''''' function density(x) as float ' atmospheric density function ' input ' x = altitude (meters) ' output ' density = atmospheric density (kilograms/cubic meter) '''''''''''''''''''''''''''''''''''''''''''''''''''''''' density = (1.0 - 0.000022556913 * x) ^ 4.256116 end function Here's a scanned PDF document for the first QuickBASIC version written in 1982. Those who remember will probably snicker at the dot matrix printout. 2017-03-19_174232_rocket1.pdf |
||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
ROCKET2 is an MMBASIC program which determines the flight performance of a single stage model rocket by numerically integrating the equations of motion. The user can specify a launch angle and the program solves the problem of non-vertical motion. The user can also display and print the results at each integration step. The information displayed includes the flight time, vertical altitude, horizontal range, velocity, mass, thrust and aerodynamic drag. This program models the variation of density with altitude and the changes in thrust and mass with time. Aerodynamic drag is updated as a function of altitude and velocity. Here's a typical user interaction with the software. program rocket2 --------------- please input the launch site altitude (meters) ? 0 please input the launch site temperature (degrees f) ? 59 please input the maximum thrust (newtons) ? 13 please input the sustainer thrust (newtons) ? 3.5 please input the time of maximum thrust (seconds) ? .14 please input the time of sustainer thrust (seconds) ? .22 please input the thrust duration (seconds) ? 1.2 please input the initial mass (grams) ? 40 please input the propellant mass (grams) ? 8.33 please input the frontal diameter (millimeters) ? 18 please input the drag coefficient (non-dimensional) ? .321 please input the launch angle (degrees [ 0 to 90 ] - measured from the ground) (for example, a vertical launch is an angle of 90 degrees) ? 90 please input the boost step size (seconds) (a value of 0.01 seconds is recommended) ? .01 please input the coast step size (seconds) (a value of 0.1 seconds is recommended) ? .1 integrating equations of motion ... program rocket2 --------------- burnout altitude 80.8596 meters burnout velocity 118.7054 meters per second coast time 8.0000 seconds total flight time 9.2000 seconds maximum altitude 459.0167 meters Here's the source code listing. ' program rocket2.bas April 2, 2017 ' determines flight performance of model rockets ' altitude at burnout (meters) ' velocity at burnout (meters per second) ' coast time and total flight time (seconds) ' maximum altitude (meters) ' numerical integration method ' MicroMite eXtreme version ''''''''''''''''''''''''''' option default float option base 1 const gravity = 9.80665 const rhosl = 1.22557 const dtr = pi / 180.0 const rtd = 180.0 / pi ' integrator coefficients const a0 = 0.045 const b0 = 0.3 const c0 = 13.0 / 126.0 const d0 = 5.0 / 18.0 const e0 = 5.0 / 42.0 const f0 = 7.0 / 600.0 const g0 = 7.0 / 30.0 const h0 = 7.0 / 15.0 const i0 = 7.0 / 6.0 const j0 = 25.0 / 63.0 const k0 = 0.7 const l0 = 19.0 / 78.0 const m0 = 35.0 / 312.0 const n0 = 15.0 / 104.0 const o0 = 64.0 / 39.0 const p0 = 70.0 / 39.0 const q0 = 15.0 / 13.0 ' global arrays and variables dim d1(2), d2(2), d3(2), d4(2) dim acc(2), xp(2), xs(2), vp(2), vs(2) dim langle, fltflg%, f5, maxthr, susthr, vex, vex2 dim massi, mprop, tt, tp, dt, tsus, tmax, tdur, t6, t7 dim thrust, mass, drag, ffactor ' begin simulation print " " print "program rocket2" print "---------------" do print " " print "please input the launch site altitude (meters)" input altsite loop until (altsite >= 0.0) do print " " print "please input the launch site temperature (degrees f)" input tempsite loop until (tempsite >= 0.0) do print " " print "please input the maximum thrust (newtons)" input maxthr loop until (maxthr > 0.0) do print " " print "please input the sustainer thrust (newtons)" input susthr loop until (susthr > 0.0) do print " " print "please input the time of maximum thrust (seconds)" input tmax loop until (tmax > 0.0) do print " " print "please input the time of sustainer thrust (seconds)" input tsus loop until (tsus > 0.0) do print " " print "please input the thrust duration (seconds)" input tdur loop until (tdur > 0.0) do print " " print "please input the initial mass (grams)" input massi loop until (massi > 0.0) massi = 0.001 * massi do print " " print "please input the propellant mass (grams)" input mprop loop until (mprop > 0.0) mprop = 0.001 * mprop do print " " print "please input the frontal diameter (millimeters)" input fdia loop until (fdia > 0.0) fdia = pi * fdia * fdia / 4000000.0 do print " " print "please input the drag coefficient (non-dimensional)" input cd loop until (cd > 0.0) do print " " print "please input the launch angle" print "(degrees [ 0 to 90 ] - measured from the ground)" print "(for example, a vertical launch is an angle of 90 degrees)" input langle loop until (langle >= 0.0 and langle <= 90.0) langle = langle * dtr do print " " print "please input the boost step size (seconds)" print "(a value of 0.01 seconds is recommended)" input dtboost loop until (dtboost > 0.0) do print " " print "please input the coast step size (seconds)" print "(a value of 0.1 seconds is recommended)" input dtcoast loop until (dtcoast > 0.0) print " " print "integrating equations of motion ..." print " " ' compensate for launch site altitude and temperature rhosite = rhosl * fna(altsite) / (1.0 + (tempsite - 59.0) / 518.67) f5 = maxthr - susthr ffactor = 0.5 * rhosite * fdia * cd t6 = tsus - tmax ' calculate exhaust velocity vex = (maxthr * 0.5 * tsus + susthr * (tdur - 0.5 * tmax - 0.5 * tsus)) / mprop vex2 = 2.0 * vex tp = 0.0 xp(1) = 0.0 xp(2) = 0.0 vp(1) = 0.0 vp(2) = 0.0 fltflg% = 0 ' powered flight np% = tdur / dtboost dt = dtboost for ii% = 1 to np% integrate tt = tp plmass plthrust ' un-comment the following line to display the ' flight conditions during the thrusting flight ' pdata next ii% ' burnout conditions xbo = xp(2) vbo = sqr(vp(1) ^ 2 + vp(2) ^ 2) ' coasting flight dt = dtcoast fltflg% = 1 do if (vp(2) > 0.0) then integrate tt = tp plmass plthrust ' un-comment the following line to display the ' flight conditions during the coasting flight ' pdata else exit do end if loop tcoast = tp - tdur ' print final conditions print " " print "program rocket2" print "---------------" print " " print "burnout altitude ", str$(xbo, 0, 4), " meters" print " " print "burnout velocity ", str$(vbo, 0, 4), " meters per second" print " " print " " print "coast time ", str$(tcoast, 0, 4), " seconds" print " " print "total flight time ", str$(tp, 0, 4), " seconds" print " " print " " print "maximum altitude ", str$(xp(2), 0, 4), " meters" print " " end '''''''''''''''''''''' '''''''''''''''''''''' function fna(x) as float ' atmospheric density function '''''''''''''''''''''''''''''' fna = (1.0 - 0.000022556913 * x) ^ 4.256116 end function '''''''' '''''''' sub diffeq ' two-dimensional equations of motion subroutine '''''''''''''''''''''''''''''''''''''''''''''''' vmag = sqr(vs(1) ^ 2 + vs(2) ^ 2) if (xp(2) < 1.0) then uvx = cos(langle) uvy = sin(langle) else uvx = vs(1) / vmag uvy = vs(2) / vmag end if drag = vmag ^ 2 * ffactor * fna(xs(2)) mass = massi - mprop thrust = 0.0 if (fltflg% = 0) then t7 = tt - tmax plmass plthrust end if tacc = (thrust - drag) / mass acc(1) = tacc * uvx acc(2) = tacc * uvy - gravity end sub ''''''''''' ''''''''''' sub integrate ' numerical integration subroutine '''''''''''''''''''''''''''''''''' tt = tp for i% = 1 to 2 xs(i%) = xp(i%) vs(i%) = vp(i%) next i% diffeq tt = tp + b0 * dt for i% = 1 to 2 d1(i%) = dt * acc(i%) xs(i%) = xp(i%) + dt * (b0 * vp(i%) + a0 * d1(i%)) vs(i%) = vp(i%) + b0 * d1(i%) next i% diffeq tt = tp + k0 * dt for i% = 1 to 2 d2(i%) = dt * acc(i%) xs(i%) = xp(i%) + dt * (k0 * vp(i%) + f0 * d1(i%) + g0 * d2(i%)) vs(i%) = vp(i%) - h0 * d1(i%) + i0 * d2(i%) next i% diffeq tt = tp + dt for i% = 1 to 2 d3(i%) = dt * acc(i%) xs(i%) = xp(i%) + dt * (vp(i%) + l0 * d1(i%) + m0 * d2(i%) + n0 * d3(i%)) vs(i%) = vp(i%) + o0 * d1(i%) - p0 * d2(i%) + q0 * d3(i%) next i% diffeq tp = tp + dt for i% = 1 to 2 d4(i%) = dt * acc(i%) xp(i%) = xp(i%) + dt * (vp(i%) + c0 * d1(i%) + d0 * d2(i%) + e0 * d3(i%)) vp(i%) = vp(i%) + c0 * (d1(i%) + d4(i%)) + j0 * (d2(i%) + d3(i%)) next i% end sub ''''''' ''''''' sub pdata ' print flight conditions subroutine '''''''''''''''''''''''''''''''''''' vmag = sqr(vp(1) ^ 2 + vp(2) ^ 2) rmag = sqr(xp(1) ^ 2 + xp(2) ^ 2) fpa = rtd * atn(vp(2) / vp(1)) accm = sgn(acc(2)) * sqr(acc(1) ^ 2 + acc(2) ^ 2) print " " print "flight time (seconds)"; tp print " " print " " print "altitude (meters)"; xp(2) print " " print "horizontal range (meters)"; xp(1) print " " print "velocity (meters per second)"; vmag print " " print "flight path angle (degrees)"; fpa print " " print "acceleration (meters/second/second)"; accm print " " print " " print "mass (grams)"; 1000.0 * mass print " " print "thrust (newtons)"; thrust print " " print "drag (newtons)"; drag print " " end sub '''''''' '''''''' sub plmass ' piecewise-linear mass subroutine '''''''''''''''''''''''''''''''''' if (tp <= tmax) then mass = massi - maxthr * tt ^ 2 / (vex2 * tmax) exit sub elseif (tp > tmax and tp <= tsus) then mass = massi - (maxthr / vex) * (tt - 0.5 * tmax) + f5 * t7 ^ 2 / (vex2 * t6) exit sub elseif (tp > tsus and tp < tdur) then mass = massi - (maxthr * tsus + susthr * t6) / vex2 - (susthr / vex) * (tt - tsus) exit sub elseif (tp >= tdur) then mass = massi - mprop exit sub end if end sub '''''''''' '''''''''' sub plthrust ' piecewise-linear thrust subroutine '''''''''''''''''''''''''''''''''''' if (tp <= tmax) then thrust = maxthr * tt / tmax exit sub elseif (tp > tmax and tp <= tsus) then thrust = maxthr - f5 * (tt - tmax) / t6 exit sub elseif (tp > tsus and tp < tdur) then thrust = susthr exit sub elseif (tp >= tdur) then thrust = 0.0 exit sub end if end sub |
||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
This MMBASIC computer program, rocket3.bas, determines the "best" initial mass of a single stage model rocket that maximizes the final altitude. It assumes a vertical launch and and compensates for non-standard launch conditions. The following PDF document provides additional information about this technique. 2017-04-03_102438_rocket.pdf Here's a typical user interaction with the software and the results calculated by the rocket3 program. program rocket3 --------------- launch site altitude (meters) ? 100 launch site temperature (degrees fahrenheit) ? 70 thrust duration (seconds) ? 1.2 total impulse (newton-seconds) ? 5 propellant mass (grams) ? 8.33 frontal diameter (millimeters) ? 18 drag coefficient (non-dimensional) ? .321 program rocket3 --------------- initial mass 24.8724 grams burnout altitude 125.1150 meters burnout velocity 190.3544 meters per second burnout mass 16.5424 grams coast time 7.5214 seconds total flight time 8.7214 seconds maximum altitude 546.2455 meters ' program rocket3.bas April 3, 2017 ' for a given model rocket engine and aerodynamic characteristics, ' this program determines the optimal launch mass of a single-stage ' model rocket which maximizes total altitude (Bengen's maxima) ' MicroMite eXtreme version ''''''''''''''''''''''''''' option default float option base 1 ' global constants const gravity = 9.80665 const rhosl = 1.22557 ' global variables dim rho, tduration, impulse, mprop, fd, cd dim altbo, velbo, massbo, weightbo dim tcoast, altcoast, tflight ' begin simulation print " " print "program rocket3" print "---------------" do print " " print "launch site altitude (meters)" input altsite loop until (altsite >= 0.0) do print " " print "launch site temperature (degrees fahrenheit)" input tempsite loop until (tempsite >= 0.0) do print "thrust duration (seconds)" input tduration loop until (tduration > 0.0) do print "total impulse (newton-seconds)" input impulse loop until (impulse > 0.0) do print "propellant mass (grams)" input mprop loop until (mprop > 0.0) do print "frontal diameter (millimeters)" input fd loop until (fd > 0.0) do print "drag coefficient (non-dimensional)" input cd loop until (cd > 0.0) ' convert mass to kilograms and diameter to square meters mprop = 0.001 * mprop fd = pi * fd ^ 2 / 4000000.0 ' compensate for launch site altitude and temperature rho = rhosl * fna(altsite) / (1.0 + (tempsite - 59.0) / 518.67) ' determine "best" initial launch mass m1 = 1.001 * mprop m2 = 10.0 * mprop tolm = 1.0e-8 minima(m1, m2, tolm, massi, altmax) ' print results print " " print "program rocket3" print "---------------" print " " print "initial mass ", str$(1000.0 * massi, 0, 4), " grams" print " " print "burnout altitude ", str$(altbo, 0, 4), " meters" print " " print "burnout velocity ", str$(velbo, 0, 4), " meters per second" print " " print "burnout mass ", str$(1000.0 * (massi - mprop), 0, 4), " grams" print " " print " " print "coast time ", str$(tcoast, 0, 4), " seconds" print " " print "total flight time ", str$(tflight, 0, 4), " seconds" print " " print " " print "maximum altitude ", str$(-altmax, 0, 4), " meters" print " " end '''''''''''''''''''''' '''''''''''''''''''''' function fna(x) as float ' atmospheric density function '''''''''''''''''''''''''''''' fna = (1.0 - 0.000022556913 * x) ^ 4.256116 end function ''''''''''''''''''''''''' ''''''''''''''''''''''''' sub alt_func(massi, altmax) ' calculate maximum altitude subroutine ''''''''''''''''''''''''''''''''''''''' thrust = impulse / tduration mass = massi - 0.5 * mprop k2 = 0.5 * rho * fd * cd weight = mass * gravity b = tduration * sqr(k2 * (thrust - weight)) / mass c = exp(b) d = exp(-b) e = 0.5 * (c + d) f = (c - d) / (c + d) ' burnout conditions altbo = (mass / k2) * log(e) velbo = f * sqr((thrust - weight) / k2) massbo = massi - mprop weightbo = massbo * gravity ' coast conditions tcoast = sqr(massbo / (k2 * gravity)) * atn(velbo * sqr(k2 / weightbo)) altcoast = (massbo / (2.0 * k2)) * log(k2 * velbo * velbo / weightbo + 1.0) ' total flight time and maximum altitude tflight = tduration + tcoast altmax = -(altbo + altcoast) end sub '''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''' sub minima(a, b, tolm, xmin, fmin) ' one-dimensional minimization ' Brent's method ' input ' a = initial x search value ' b = final x search value ' tolm = convergence criterion ' output ' xmin = minimum x value ' note ' user-defined objective subroutine ' coded as alt_func(x, fx) ' remember: a maximum is simply a minimum ' with a negative attitude! ''''''''''''''''''''''''''''''''''''' ' machine epsilon LOCAL epsm = 2.23e-16 ' golden number LOCAL c = 0.38196601125 LOCAL d, e, t2, p, q local r, u, fu LOCAL x, xm, w local v, fx, fw LOCAL tol1, fv x = a + c * (b - a) w = x v = w e = 0.0 p = 0.0 q = 0.0 r = 0.0 alt_func(x, fx) fw = fx fv = fw for iter% = 1 to 100 if (iter% > 50) then print ("error in subroutine minima!") print ("(more than 50 iterations)") end if xm = 0.5 * (a + b) tol1 = tolm * abs(x) + epsm t2 = 2.0 * tol1 if (abs(x - xm) <= (t2 - 0.5 * (b - a))) then xmin = x fmin = fx exit sub end if if (abs(e) > tol1) then r = (x - w) * (fx - fv) q = (x - v) * (fx - fw) p = (x - v) * q - (x - w) * r q = 2.0 * (q - r) if (q > 0.0) then p = -p end if q = abs(q) r = e e = d end if if ((abs(p) >= abs(0.5 * q * r)) or (p <= q * (a - x)) or (p >= q * (b - x))) then if (x >= xm) then e = a - x else e = b - x end if d = c * e else d = p / q u = x + d if ((u - a) < t2) or ((b - u) < t2) then d = sgn(xm - x) * tol1 end if end if if (abs(d) >= tol1) then u = x + d else u = x + sgn(d) * tol1 end if alt_func(u, fu) if (fu <= fx) then if (u >= x) then a = x else b = x end if v = w fv = fw w = x fw = fx x = u fx = fu else if (u < x) then a = u else b = u end if if ((fu <= fw) or (w = x)) then v = w fv = fw w = u fw = fu elseif ((fu <= fv) Or (v = x) Or (v = w)) then v = u fv = fu end if end if next iter% end sub |
||||
Gizmo Admin Group Joined: 05/06/2004 Location: AustraliaPosts: 5026 |
Thanks cdeagle, nice work. Speaking of rockets and software, it would be fascinating to talk to the software engineers who work on the SpaceX Falcon 9 first stage. That thing is a engineering masterpiece. Glenn The best time to plant a tree was twenty years ago, the second best time is right now. JAQ |
||||
Print this page |