[Python] Simulating rocket trajectory to ISS

[codymc94]

background information is that I have been working on code for a rocket launch to the ISS and I have gotten it close. The problem is when calculating the Force of drag, the problem occurs with it. As entering space, and for this "model" were going to say that temperature is 0 K in space or basically just out of the atmosphere. Obviously absolute 0 isn't possible, but perhaps that would be a simple solution to the problem. When entering space, T= 0 K, and there isn't any air, well not in the "model". So since there isn't any air then there wouldn't be any drag force. I did the calculation for my formula and at y= 44307.69 m, the rocket leaves the atmosphere. the equations I use for temperature is T=Tz - L*y where Tz is the Temperature at time 0, L is temperature lapse rate (assuming linear as y increases) of value 0.0065 K/m, and y is height in m.

Next, is the pressure, and as you know as y increases pressure should also go down. Now as leaving the atmosphere there should not be any pressure. This formula is P= Pnot * (T/Tz)^((g*M)/(R*L)); where Pnot is pressure @ y=0 or 101 kPa, T=current temperature, Tz= Temp @ y=0 (288K), g=gravity (which uses another formula but I have that one working correctly), M= molar mass of air value of 0.02896 kg/Mole, R=8.314 J/K*M, and L is the temperature lapse rate again.

Anyways this is a lot of information obviously, and my question is how to get the force of drag to stop at that y= 44307.69 m, and then keep calculating in the while loop for the computer program. also using the Euler method, in the python language.

Trying to get it to go to the ISS and then come back to earth.

Any help would be appreciated. Thanks.COPY OF THE CODE ATTACHED BELOW

from __future__ import division, print_function
from pylab import *
from scipy.interpolate import interp1d

dt= 0.01 #seconds, the timestep
t= 0.0 #seconds, initial time

yfinal= 355000.0 #m, =355km how far up the ISS is traveling
y= 0.0 #m at ground level
v= 0.0 #m/s initial speed
pay= 5000.0 #kg
me= 10000.0 #kg, empty mass with NO fuel or payload
gas= 10000.0 #kg of rocket fuel, 10,000 is max fuel
m= 20000.0 #kg   ###+ pay #total mass, empty rocket, fuel, and payload
Em= 10000.0
Ev= 7000.0 #m/s exhaust velocity
Br= 100.0 #kg/s  burn rate of the fuel
A= 10.0 #m**2   Reference or cross sectional area
Cd= 0.5  #drag coefficient

gs=9.8 #m/s**2  initial gravity at the Earth's surface
Vex= 7000.0 #m/s exit velocity of the fuel being ejected out
Ft= Br * Vex #force of thrust, acceleration upward
a=0.0
nogas = gas/Br  #the point at which the rocket runs out of gas
Re= 6378100.0 #meters radius of the earth
g= 0.0
rho= 0.0
R= 8.314 #Joules/Kelvin*Mole
M= 0.02896 #kg/mole        
T= 288.0
Pz=101000.0  #Pascal  Pressure at time 0
Tz= 288.0 #Kelvin   Tempature at time 0. Aprox. 14.85 Celsius
L= 0.0065 #Kelvin/meter  "Tempature lapse rate"
P= 101000.0 #Pressure in Pascal
Fd= 0.0
Fg= 196000
#initial array values
tArray = array([t])
yArray = array([y])
vArray = array([v])
aArray= array([a]) 
gArray= array([g])  
mArray= array([m])
FtArray= array([Ft])
FdArray = array([Fd])
FgArray= array([Fg])    
TArray = array([T])   
rhoArray= array([rho])
PArray = array([P])
RArray= array([R])
MArray= array([M])
#calculations for arrays
#if y<=44307.69231: 
   # break(T)
while (y<=yfinal):# and (T>=0):# and (P>=10100): # and (nogas>=t):
    
    t= t + dt  #increments TIME
    
    gas= gas - Br*dt
    m = Em + gas
    Ft= Br * Vex

    T= Tz - (L*y)#*dt
    
    g=  (gs/((1+y/Re)**2))
    P= Pz*((abs(T/Tz))**((g*M)/(R*L)))
    rho= (P*M)/(R*T) #*dt
    if(T < 0):
        Fd= 0 
    
    Fd=  .5*rho*A*Cd*v*abs(v)
    
    
    Fg= (m)*g
    a=(-Fg-Fd+Ft)/m
    
    y= y + v*dt #increments HEIGHT
         
    v= v + a*dt #increments VELOCITY

       
    
    yArray = append(yArray, y) #fills in height array
    vArray = append(vArray, v)#fills in velocity array
    tArray = append(tArray, t)#fills in time array
    aArray = append(aArray, a)
    gArray = append(gArray, g)
    mArray = append(mArray, m)
    FtArray= append(FtArray, Ft)
    FgArray= append(FgArray, Fg)
    FdArray= append(FdArray, Fd)
    TArray = append(TArray, T)
    rhoArray= append(rhoArray, rho)
    PArray = append(PArray, P)
    RArray = append(RArray, R)
    MArray= append(MArray, M)
    print("time", t, "secs", '\t', "height", y, "m", '\t', 'velocity', v, "m/s")
    print(aArray)

 

[mfb]

Anyways this is a lot of information obviously, and my question is how to get the force of drag to stop at that y= 44307.69 m, and then keep calculating in the while loop for the computer program. also using the Euler method, in the python language.

Your approach with if is good:

if(condition)
   Fd= 0 
else
   Fd= .5*rho*A*Cd*v*abs(v)

There are some other issues, however:

L is temperature lapse rate (assuming linear as y increases) of value 0.0065 K/m

This is true close to the ground, maybe for the first 10km, but it is completely wrong afterwards. Your atmospheric model is not useful - above the height where your atmosphere ends, the space shuttle performed its re-entry maneuvers, and it was glowing red from the drag there.

Your model seems to be one-dimensional, looking only at the height. That is completely unrealistic - the ISS is orbiting Earth with a speed of ~7-8km/s, and most of the rocket acceleration is needed to get this speed. The height is not so crucial, it is just needed to get out of the atmosphere.

 

[D H]

The ISS is *not* out of the atmosphere. Its orbit decays due to atmospheric drag. Here's a plot of the ISS's altitude from 1999 to 2009.

The discontinuities are where the ISS did a reboost to raise its altitude. The slopes of the continuous segments vary quite a bit. Notice how the ISS maintained a rather high altitude during 2001-2002. The periods between reboosts are rather steep lines. Notice how the altitude is much lower in 2006-2007, and yet even with that significant decrease in altitude, the drag losses in 2006-2007 are much lower than they were in 2001-2002.

I wouldn't expect you to model this, but the reason is the solar cycle. Solar max puffs up the upper atmosphere. Solar cycle 23 had a maximum in 2000 and another in 2002. The ISS altitude was raised during that solar max so as to reduce the need for reboosts. (That sharp drop in altitude in April-May 2000: That's because of STS-101, whose launch was delayed a number of times. NASA didn't want to do a reboost during that interval unless they absolutely had to.)

What about the flip side, the period encompassing 2006 and 2007? Launches are considerably cheaper if the launch vehicle doesn't have to go so high. Solar cycle 23 ended with a prolonged minimum, and the Earth's upper atmosphere contracted in response. NASA and Russia took advantage of this extended lull by keeping the ISS altitude low during 2006 and 2007.

While NASA and Russia do need to model this behavior, you don't. That takes a fairly sophisticated atmosphere model. Having an atmosphere model that has the atmosphere just end is not a great idea. Having an atmosphere model with a constant lapse rate until the atmosphere ends is a very bad idea. Going from sea level on up, temperature in the Earth's atmosphere drops with increasing altitude in the troposphere, then rises in the stratosphere, then drops again in the mesosphere, and rises once again in the thermosphere. A fairly simple model that exhibits these behaviors is the US Standard Atmosphere model. There's even a python model that incorporates the US Standard Atmosphere, up to 84 km altitude: https://pypi.python.org/pypi/AeroCalc.