Newer
Older
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
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
if T == T0 : # isohermal region
P = P0*np.exp(-g0*(altitude-h0_tab[index])/(r*T))
else : # gradient region
P = P0*np.power(T/T0,-g0/(a*r))
rho = P/(r*T)
return T,P,rho
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)
if Aadapted > Ae: # the nozzle is under-expanded
# compute the exit Mach number iteratively using dichotomy algorithm
A_ratio_target = Ae/Ath
Me_upper = Madapted
Me_lower = 1
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_upper = 1
Me_lower = 0
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)
Pe = Pt*np.power(Tt/Te,gamma/(1-gamma))
while (Pe-Pamb)*1e-5 > tolerance:
Mshock_old = Mshock
Mshock += step
Pe_old = Pe
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_upper = 1
Me_lower = 0
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:
Loading full blame...