Newer
Older
# -*- coding: utf-8 -*-
"""
Created on Sun May 12 20:08:52 2019
@author: MAC
"""
import numpy as np
from numpy import linalg as LA
def initial_vel(latitude,retrograde=False):
"""
Initial velocity (m.s-1) depending on the launch point latitude (°)
"""
T_side = 86164.10 # sidereal day (s)
R_T = 6378e3 # earth raduis (m)
initial_velocity = (2*np.pi/T_side)*R_T*np.cos(np.deg2rad(np.abs(latitude)))
if retrograde:
return -initial_velocity
else:
return initial_velocity
def final_target_vel(h, a):
"""
inputs:
h: injection altitude (m)
a: orbit semi major axis (m)
output:
orbit absolute injection velocity (m.s-1)
"""
mu = 3.9860e14 # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2)
R_T = 6378e3 # earth raduis (m)
return np.sqrt(2*mu*((1/(R_T+h))-(1/(2*a))))
def staging_optim(isp1,isp2,eps1,eps2,Vc,Mp):
"""
Optimise mass loading in the case of two stage rocket
inputs:
isp1, isp2 : specific impulse of the 1st and 2nd stage (s)
eps1, eps2 : structural coefficient of the 1st and 2nd stage (adimentional)
Vc : target critical velocity (m.s-1)
Mp : payload mass (Kg)
outputs:
M1 : mass of the 1st stage
M2 : mass of the 2nd stage
NB:
structure mass
structural coefficient = --------------------------------
structure mass + propellant mass
ctritical velocity = final velocity - initial velocity + losses
"""
# Notation :
#
# Mt = total mass
# q1 = M1/Mt
# q2 = M2/Mt
# qp = Mp/Mt
# omega = qp + q2
g0 = 9.81 # standard gravity (m.s-2)
omega_tab=np.linspace(0.001,0.9,100) # array of omega to test
qp_tab=[]
for omega in omega_tab:
q2 = omega+eps1*(1-omega)
q2 = np.power(q2,isp1/isp2)*np.exp(Vc/(g0*isp2))
q2 = 1/q2
q2 = 1-q2
q2 *= omega/(1-eps2)
qp = omega-q2
qp_tab.append(qp)
qp=max(qp_tab)
omega=omega_tab[qp_tab.index(qp)]
q2=omega-qp
q1=1-q2-qp
M1=(Mp/qp)*q1
M2=(Mp/qp)*q2
return M1, M2
def standard_atmosphere(altitude):
"""
Compute atmospheric properties according to the International Standaard Atmosphere (ISA) model
input:
altitude (m)
outputs:
temperature (K)
pressure (Pa)
density (Kg.m-3)
assert altitude < 84852 , str(altitude)+' is higher than max altitude of the ISA model'
g0 = 9.81 # standard gravity (m.s-2)
r = 287 # specific gas constant (J.Kg-1.K-1)
# ISA model
h0_tab = np.array([0, 11e3, 20e3, 32e3, 47e3, 51e3, 71e3, 84852]) # base altitude (m)
a_tab = np.array([-6.5e-3, 0, 1e-3, 2.8e-3, 0, -2.8e-3, -2.0e-3]) # lapse rate (K.m-1)
# the base properties of each range are tabulated to reduce computational cost
T0_tab = np.array([288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65]) # base temperature (K)
P0_tab = np.array([101325, 22632, 5474.9, 868.02, 110.91, 66.939, 3.9564]) # base pressure (Pa)
index = 0
while altitude > h0_tab[index+1]:
index += 1
a = a_tab[index]
T0 = T0_tab[index]
P0 = P0_tab[index]
T = T0+a*(altitude-h0_tab[index]) # Toussaint equation
P = P0*np.exp(-g0*(altitude-h0_tab[index])/(r*T))
P = P0*np.power(T/T0,-g0/(a*r))
rho = P/(r*T)
def laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,r,tolerance=1e-5):
"""
Compute the thrust assuming a started nozzle
inputs:
Pt : total pressure in the combustion chambre (Pa)
Tt : total temperature in the combustion chambre (K)
Pamb : ambient pressure (Pa)
Ath : throat section area (m2)
Ae : exit section area (m2)
gamma : ratio of specific heats (adimentional)
r : specific gas constant (J.Kg-1.K-1)
tolerance :
outputs:
Thrust (N)
Specific impulse (s)
N.B : Thrust = mass flow rate * exit velocity + (exit pressure - ambiant pressure) * exit area
"""
assert gamma > 1, 'gamma must be > 1'
g0 = 9.81 # standard gravity (m.s-2)
A_ratio = lambda M,gamma: (1/M)*np.power((2/(gamma+1))*(1+((gamma-1)/2)*M*M),(gamma+1)/(2*gamma-2))
# compute the nozzle exit section area corresponding to the adapted regime
Madapted = np.sqrt((2/(gamma-1))*(np.power(Pamb/Pt,(1-gamma)/gamma)-1))
Aadapted = Ath*A_ratio(Madapted,gamma)
print(Aadapted,Ae)
if Aadapted > Ae: # the nozzle is under-expanded
# compute the exit Mach number iteratively using dichotomy algorithm
A_ratio_target = Ae/Ath
Me = (Me_upper+Me_lower)/2
error = A_ratio_target-A_ratio(Me,gamma)
while np.abs(error) > tolerance:
if error > 0:
Me_lower = Me
else:
Me_upper = Me
Me = (Me_upper+Me_lower)/2
error = A_ratio_target-A_ratio(Me,gamma)
Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature
Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure
elif Aadapted < Ae: # the nozzle is over-expanded
Mshock = 1
step = 0.1
A_ratio_target = Ae/Ath
Me = (Me_upper+Me_lower)/2
error = A_ratio_target-A_ratio(Me,gamma)
while np.abs(error) > tolerance:
if error < 0:
Me_lower = Me
else:
Me_upper = Me
Me = (Me_upper+Me_lower)/2
error = A_ratio_target-A_ratio(Me,gamma)
print("erreur: ",error)
print("Mach", Me)
Te = Tt/(1+((gamma-1)/2)*Me*Me)
Pe = Pt*np.power(Tt/Te,gamma/(1-gamma))
while (Pe-Pamb)*1e-5 > tolerance:
Mshock += step
Pt_ratio = np.power(((gamma+1)*Mshock*Mshock)/((gamma-1)*Mshock*Mshock+2),gamma/(gamma-1))*np.power((gamma+1)/(2*gamma*Mshock*Mshock-gamma+1),1/(gamma-1))
Pi = Pt*Pt_ratio
Ai = Ath/Pt_ratio
A_ratio_target = Ae/Ai
Me = (Me_upper+Me_lower)/2
error = A_ratio_target-A_ratio(Me,gamma)
while np.abs(error) > tolerance:
if error < 0:
Me_lower = Me
else:
Me_upper = Me
Me = (Me_upper+Me_lower)/2
error = A_ratio_target-A_ratio(Me,gamma)
#print("erreur 1: ",error)
#print("Mach 1:", Me)
Te = Tt/(1+((gamma-1)/2)*Me*Me)
Pe = Pi*np.power(Tt/Te,gamma/(1-gamma))
if (Pe > Pamb) and (A_ratio(Mshock,gamma) > Ae/Ath):
# the normal shock is outside the divergent and all the flow inside is supersonic
A_ratio_target = Ae/Ath
Me = (Me_upper+Me_lower)/2
error = A_ratio_target-A_ratio(Me,gamma)
while np.abs(error) > tolerance:
if error > 0:
Me_lower = Me
else:
Me_upper = Me
Me = (Me_upper+Me_lower)/2
error = A_ratio_target-A_ratio(Me,gamma)
Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature
Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure
break
if Pe < Pamb:
Pe = Pe_old
Mshock = Mshock_old
step /= 2
else: # the nozzle is adapted
Me = Madapted
Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature
Ve = Me*np.sqrt(gamma*r*Te) # exit velocity
# maximum mass flow rate in the case of a started nozzle
Q = np.power(2/(gamma+1),(gamma+1)/(2*(gamma-1)))*np.sqrt(gamma/(Tt*r))*Pt*Ath
Thrust = Q*Ve+(Pe-Pamb)*Ae
isp = Thrust/(Q*g0)
return Thrust,isp
def adiabatic_flame(mech,T,P,X):
"""
Compute adiatabic flame temperature
inputs:
mech : chemical mechanism in cti format
T : initial temperature (K)
P : pressure (Pa)
X : initial chemical composition
def injector():
"""
"""
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
#%% Trajectory tools
def solve_traj():
"""
"""
# .......
#%% Astrodynamics tools
def state_to_orbital_elem(r,v):
"""
Convert state vector to orbital elements
inputs:
r : coordinate vector (m)
v : velocity vector (m.s-1)
output: {'a','e','i','capital_omega','small_omega','nu'}
where:
a : semi_major axis (m)
e : eccentricity (adimentional)
i : inclination (deg)
capital_omega : argument of the ascending node (deg)
small_omega : argument of the perigee (deg)
nu : true anomaly (deg)
NB: - An internal quantity (small) is considered for the case of undefined
classical orbital elements. When an element is undefined it is set to
zero and other(s) following it are measured according to this convention
- The function stops when r and v are aligned
- This function is widely inspired by the Matlab function written by
Prof. Josep Masdemart from UPC
"""
orbel={'a':0,'e':0,'i':0,'capital_omega':0,'small_omega':0,'nu':0}
mu = 3.9860e14 # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2)
small = 1e-10
h = np.cross(r,v) # angular momentum
rm = LA.norm(r) # module of the coordinate vector
vm = LA.norm(v) # module of the velocity vector
hm = LA.norm(h) # module of the angular momentum vector
if hm < small:
print('rectilinar or collision trajectory')
return
n = np.array([-h[1],h[0],0])
nm = LA.norm(n)
if nm < small:
n = np.array([1,0,0])
nm = 0
else:
n = n/nm
e = ((vm*vm-mu/rm)*r-dotrv*v)/mu
em = LA.norm(e)
if em < small:
e = n
em = 0
else:
e = e/em
orbel['a'] = (hm*hm/mu)/(1-em*em)
aux = h[2]/hm
if np.abs(aux) > 1: aux = np.sign(aux)
# all the angles are first computed in rad
orbel['i'] = np.arccos(aux)
if nm > small:
aux = n[0]
if np.abs(aux) > 1: aux = np.sign(aux)
orbel['capital_omega'] = np.arccos(aux)
if n[1] < 0: orbel['capital_omega'] = 2*np.pi-orbel['capital_omega']
if em > small:
aux = np.dot(e,n)
if np.abs(aux) > 1: aux = np.sign(aux)
orbel['small_omega'] = np.arccos(aux)
if e[2] < 0: orbel['small_omega'] = 2*np.pi-orbel['small_omega']
aux = np.dot(e,r)/rm
if np.abs(aux) > 1: aux = np.sign(aux)
orbel['nu'] = np.arccos(aux)
if em > small:
if dotrv < 0: orbel['nu'] = 2*np.pi-orbel['nu']
if np.dot(h,ha) < 0: orbel['nu'] = 2*np.pi-orbel['nu']
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
# convert the angles to deg
orbel['i'] = np.rad2deg(orbel['i'])
orbel['capital_omega'] = np.rad2deg(orbel['capital_omega'])
orbel['small_omega'] = np.rad2deg(orbel['small_omega'])
orbel['nu'] = np.rad2deg(orbel['nu'])
return orbel
def propagate_orbit():
"""
"""
#%% Thermodynamic cycles
# the following functions provide the pressure, the mixuture composition
# and the initial temperature in the combustion chambre and properties of some
# components of the cycle.
def pressure_fed():
"""
"""
def electric_turbo():
"""
"""
def closed_expander():
"""
"""
def dual_expander():
"""
"""
def open_expander():
"""
"""