Robert Heinlein
Aleksandra Janczewska,
Katarzyna Turbańska
Od czasu wynalezienia prochu strzelniczego w Chinach ponad siedem wieków temu, ludzie nauczyli się wysyłać w niebo obiekty przy pomocy kontrolowanych eksplozji. Mechanizmy zachodzące w tych statkach i ich silnikach, zwanych rakietami, wykorzystywane są przy produkcji m.in. fajerwerków, flar sygnałowych czy broni.
Ale od lat pięćdziesiątych XX wieku rakiety pozwalają nam także wysyłać satelity, sondy, roboty, zwierzęta i ludzi na orbitę wokół Ziemi - a nawet poza nią.$^{[1]}$
Chociaż eksploracja kosmosu jest oczywiście głównym celem wszystkich wysiłków naukowców, warto pamiętać, że „wyjście poza” Ziemię daje nam lepsze zrozumienie naszej własnej planety: prognozowanie pogody, badanie klimatu i nawigacja to tylko trzy aspekty, które udoskonalamy dzięki rozwojowi rakiet kosmicznych.
Cytując klasyka amerykańskiej fantastyki naukowej, Roberta Heinleina: „Jeśli potrafisz wynieść swój statek na orbitę, jesteś w połowie drogi dokądkolwiek”$^{[2]}$. Oznacza to mniej więcej tyle, że w przypadku każdej misji kosmicznej opuszczenie powierzchni Ziemi i następujące po nim swobodne poruszanie się po jej orbicie pochłania połowę kosztów energii.
Jeśli chcemy zrozumieć działanie rakiet, musimy zacząć od pojęcia kosmosu.
Lecąc rakietą w kierunku gwiazd, nie spotkamy żadnego znaku drogowego, głoszącego: „Kosmos: Populacja 0, jedź ostrożnie”. Grawitacja sięga nieskończoności, dlatego nie ma wyraźnej linii podziału między końcem Ziemi a początkiem kosmosu.
Umowna granica między atmosferą ziemską oraz przestrzenią kosmiczną znajduje się na wysokości 100 km nad Ziemią i nazywana jest Linią Kármána. Nie oznacza to wcale, że atmosfera w tym miejscu magicznie znika! W końcu satelity ziemskie latają na wysokościach dochodzących nawet 160 km, czyli ponad 10 razy wyżej niż samoloty.
Można by pomyśleć, że kosmos jest daleko, jednakże sto kilometrów to wcale nie tak dużo. Zwykłemu samochodowi przebycie takiej drogi zajęłoby godzinę; rakieta natomiast potrafi pokonać tę trasę około 20 razy szybciej - w zaledwie 3 minuty.
Z punktu widzenia kogoś, kto projektuje rakietę, przestrzeń kosmiczna jest faktycznie miejscem poza zasięgiem Ziemi, gdzie nie działa na nas grawitacja ziemska i nie ma atmosfery. Chociaż myślimy o kosmosie jako o próżni, nie jest on całkowicie pusty. Znajdziemy w nim przecież inne planety, meteoryty czy kosmiczne śmieci - pozostałości po zniszczonych rakietach bądź satelitach.
Rakieta kosmiczna jest pojazdem z silnikiem odrzutowym o ogromnej mocy, przeznaczonym do przewożenia ludzi lub sprzętu poza Ziemię. Jeśli zdefiniujemy kosmos jako przesteń poza atmosferą ziemską, to oznacza to, że nie ma tam wystarczającej ilości tlenu, by zasilić tradycyjny silnik spalinowy. Rakieta potrzebuje więc silnika, który posiada własne źródło tlenu.$^{[3]}$
Ruch rakiety oparty jest na trzeciej zasadzie dynamiki Newtona, potocznie nazywanej „prawem akcji i reakcji”. Rakieta spala paliwo i wyrzuca gorący gaz poprzez otwór wylotowy, wytwarzając w ten sposób siłę akcji, która nadaje gazowi prędkość o zwrocie przeciwnym do zwrotu prędkości rakiety. Siłę reakcji, działającej na rakietę, nazywamy ciągiem. Ciąg zwiększa prędkość rakiety, co skutkuje pojawieniem się przyspieszenia.
Rakieta - obiekt latający wprawiany w ruch silnikiem odrzutowym, służący do lotów kosmicznych.$^{[4]}$
Jeżeli ruch rakiety odbywa się w przestrzeni kosmicznej, to siły zewnętrzne ($F_{zew}$) są zaniedbywane i wtedy zmiana pędu rakiety jest równa sile ciągu (jest spełniona zasada zachowania pędu). Natomiast gdy ruch odbywa się w pobliżu Ziemi (np. tuż po starcie), to wówczas $F_{zew}$ reprezentuje ciężar rakiety oraz siłę oporu atmosfery i trzeba je uwzględnić. Konstruktorzy rakiet starają się uzyskać jak największą siłę ciągu, aby przezwyciężyć $F_{zew}$.$^{[5]}$
Rakieta wielostopniowa, rakieta wieloczłonowa - rakieta podzielona na kilka członów (stopni) napędzających ją kolejno. Każdy człon ma własny silnik i materiały pędne. Ostatni stopień zawiera statek kosmiczny lub (w przypadku pocisku rakietowego) głowicę z ładunkiem bojowym.$^{[6]}$
Zalety rakiet wielostopniowych$^{[7]}$:
Wady rakiet wielostopniowych$^{[7]}$:
Rozważmy rakietę, na którą działa siła grawiatcji, skierowana pionowo w dół. Od momentu startu ($t = 0$) rakieta wyrzuca gaz przy stałym strumieniu masy $R = \frac{dm}{dt}$ (gdzie $dm$ - masa wyrzucanego gazu, a $dt$ - czas, w którym gaz jest wyrzucany). Prędkość spalin względem rakiety jest równa $v_{e}$.
Przyspiesznie rakiety wynosi wtedy: $$a = \frac{v_{e}}{m}R - g = \frac{v_{e}}{m}\frac{dm}{dt} -g.$$
Poszczególne elementy we wzorze mają następujące znaczenie:
Prędkość końcową rakiety, a także prędkość w każdym momencie jej lotu, możemy opisać za pomocą równania Ciołkowskiego.
Równanie rakiety Ciołkowskiego – równanie matematyczne, opisujące ruchy pojazdów, zachowujących się tak jak rakiety. Są to urządzenia, które mogą nadawać sobie przyspieszenie za pomocą siły ciągu poprzez pozbycie się części swojej masy z dużą prękością i w ten sposób zachowując pęd.
\begin{equation} \Delta v = v_{e} \ln{\frac{m_{0}}{m_{1}}}, \end{equation}gdzie:
Obliczmy stosunek mas potrzebny do ucieczki spod oddziaływania grawitacji ziemskiej, zaczynając od momentu, kiedy rakieta spoczywa i zakładając, że prędkość ucieczki z Ziemi wynosi około $11,2 \cdot 10^3\frac{m}{s}$, a prędkość wyrzucania spalin względem rakiety to $v_{e}= 2,5 \cdot 10^3\frac{m}{s}$.
Przekształcamy wzór Ciołkowskiego i otrzymujemy: $$\frac{m_{0}}{m_{1}} = e^{\frac{v}{v_{e}}}.$$
Następnie podstawiamy dane: $$\frac{m_{0}}{m_{1}} = e^{\frac{11,2 \cdot 10^3\frac{m}{s}}{2,5 \cdot 10^3\frac{m}{s}}} = e^{4,48} = 88,$$ $$m_{1} = \frac{m_{0}}{88}.$$
Wynika z tego, że po spaleniu paliwa pozostaje jedynie $\frac{1}{88}$ masy początkowej, czyli paliwo stanowi, aż $\frac{87}{88}$ masy początkowej rakiety.
Rysunek poglądowy
Wiemy, że prędkość spalin względem rakiety ($v_{e}$) dana jest zależnością: $$v_{e} = v - v_{s},$$ gdzie:
Jeżeli w przedziale czasu $dt$ z rakiety zostaje wyrzucona masa $dm_{s}$ z prędkością $v_{s}$ to masa rakiety maleje o $dm$, a jej prędkość rośnie o $\Delta v$. Przy czym: $$\frac{dm_{s}}{dt} = -\frac{dm}{dt},$$ ze znakiem minus, ponieważ rakieta traci swoją masę.
Zmiana pędu układu ($dP$) w czasie $dt$ wynosi: \begin{equation} \frac{dP}{dt} = \frac{dP_{rakiety}}{dt} + \frac{dP_{spalin}}{dt} = \frac{d(mv)}{dt} + v_{s}\frac{dm_{s}}{dt} = m\frac{dv}{dt} + v\frac{dm}{dt} + v_{s}\frac{dm_{s}}{dt}, \end{equation} zakładając, że prędkość wyrzucania spalin jest stała.
Wiemy, że: $$-v_{e}\frac{dm_{s}}{dt} = v\frac{dm}{dt} + v_{s}\frac{dm_{s}}{dt}.$$
Zgodnie z drugą zasadą dynamiki Newtona: \begin{equation} F_{zew} = \frac{dP}{dt} = m\frac{dv}{dt} - v_{e}\frac{dm_{s}}{dt}. \end{equation}
Ostatni wyraz w równaniu może być interpretowany jako siła wywierana na układ przez substancję (spaliny), która z niego wylatuje. W przypadku rakiety (samolotu) nosi ona nazwę siły ciągu.
Jeżeli na układ nie działają żadne zewnętrzne siły, wtedy $F_{zew} = 0$.
Wynika z tego, że: $$m\frac{dv}{dt} = v_{e}\frac{dm_{s}}{dt}.$$
Po podstawieniu $dm_{s} = -dm$ i przekształceniu równania otrzymujemy: $$dv = -v_{e}\frac{dm}{m}.$$
Zakładając, że $v_{e}$ jest stała otrzymujemy: \begin{equation} \int_{0}^{\Delta v} \,dv = -v_{e}\int_{m_{0}}^{m_{1}} \frac{1}{m} \, dm, \end{equation} gdzie:
Ostatecznie otrzymujemy: \begin{equation} \Delta v = v_{e}\ln{\frac{m_{0}}{m_{1}}}. \end{equation}
Wzór na zmianę prędkości rakiety w czasie: $$v = v_e\ln\left( \frac{m}{m - \frac{dm}{dt}t} \right),$$ gdzie:
Drogę pokonaną przez rakietę obliczamy ze wzoru $s = v\cdot t$, podstawiając $v$ obliczone w każdej sekundzie jej lotu, więc $t$ zawsze równa się $1$s.
Przyspieszenie rakiety obliczamy ze wzoru $a = \frac{\Delta v}{\Delta t}$ dla $\Delta t = 1$s.
from IPython.display import HTML, display
def display_table(data):
"""Create html table from data table given as argument."""
html = "<table>"
for row in data:
html += "<tr>"
for field in row:
html += "<td><h4>%s</h4><td>"%(field)
html += "</tr>"
html += "</table>"
display(HTML(html))
if __name__ == "__main__":
#zmienna masa
M1 = [2500000, 2916500, 3100000, 3256500]
ve1 = [34000, 34000, 34000, 34000]
f1 = [1000, 1000, 1000, 1000]
#zmienna prędkość spalin
M2 = [2916500, 2916500, 2916500, 2916500]
ve2 = [31000, 32000, 33000, 34000]
f2 = [1000, 1000, 1000, 1000]
#zmienne tempo spalania
M3 = [2916500, 2916500, 2916500, 2916500]
ve3 = [34000, 34000, 34000, 34000]
f3 = [1000, 1200, 1400, 1600]
for j in range(1,4):
data = [["TABELA {}".format(j),"masa rakiety [kg]", "prędkość spalin [m/s]", "tempo spalania [kg/s]"]]
for i in range(len(M1)):
data.append(["rakieta {}".format(i+1), eval("M{}".format(j))[i], eval("ve{}".format(j))[i], eval("f{}".format(j))[i]])
display_table(data)
print("---------------------------------------------------------------------")
import math
import matplotlib.pyplot as plt
import numpy as np
def distance(v):
"""Count distance travelled by rocket in every second.
Args:
v [list]: rocket velocity in every second.
Return:
[list]: distance travelled by rocket in every second.
"""
s = [0]
current_s = 0
for current_velocity in v[0:len(v)-1]:
current_s += current_velocity
s.append(current_s/1000)
return s
def velocity(i, M, ve, f):
"""Count rocket velocity.
Args:
i [int]: time from launching rocket (in seconds),
M [int]: mass of rocket with propellant,
ve [int]: exhaust velocity,
f [int]: mass of ejected gas per second.
Return:
[int]: rocket velocity in i.
"""
v = ve * math.log(M/(M-f*i))
return v
def acceleration(v):
"""Count rocket acceleration.
Args:
v [list]: rocket velocity in every second.
Return:
[list]: rocket acceleration in every second.
"""
a = []
for i in range(len(v)-1):
a.append((v[i+1]-v[i])/1)
return a
def plot1(x, M, ve, f, variable, var_values, number):
"""Draw velocity vs. time plot, acceleration vs. time plot and distance vs. time plot.
Args:
x [numpy.linspace]: time interval (in seconds),
M [list]: masses of rocket with propellant,
ve [list]: exhaust velocities,
f [list]: masses of ejected gas per second.
"""
fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, sharex=True, figsize=(19, 5.5))
v0 = [velocity(i, M[0], ve[0], f[0]) for i in x]
v1 = [velocity(i, M[1], ve[1], f[1]) for i in x]
v2 = [velocity(i, M[2], ve[2], f[2]) for i in x]
v3 = [velocity(i, M[3], ve[3], f[3]) for i in x]
ax0.set_title('{}.1.Zależność drogi pokonanej przez rakietę od czasu lotu'.format(number), fontsize=13)
for list_of_velocities in [v0, v1, v2, v3]:
ax0.plot(x, distance(list_of_velocities))
ax0.legend([variable + ": " + str(v) for v in var_values])
ax0.set_ylim(0)
ax0.set_xlabel('czas [s]', fontsize=10)
ax0.set_ylabel('droga [km]', fontsize=10)
ax1.set_title('{}.2.Zależność prędkości rakiety od czasu lotu'.format(number), fontsize=13)
for list_of_velocities in [v0,v1,v2,v3]:
ax1.plot(x, list_of_velocities)
ax1.legend([variable + ": " + str(v) for v in var_values])
ax1.set_ylim(0)
ax1.set_xlabel('czas [s]', fontsize=10)
ax1.set_ylabel('prędkość [m/s]', fontsize=10)
ax2.set_title('{}.3.Zależność przyspieszenia rakiety od czasu lotu'.format(number), fontsize=13)
for list_of_velocities in [v0,v1,v2,v3]:
ax2.plot(x[1:len(x)], acceleration(list_of_velocities))
ax1.legend([variable + ": " + str(v) for v in var_values])
ax2.legend([variable + ": " + str(v) for v in var_values])
ax2.set_ylim(0)
ax2.set_xlabel('czas [s]', fontsize=10)
ax2.set_ylabel('przyspieszenie [m/s²]', fontsize=10)
fig.suptitle('Wykresy zależności, zmienna - {}'.format(variable), fontsize=16)
plt.show()
if __name__ == "__main__":
x = np.linspace(0, 1200, 1201)
plot1(x, M1, ve1, f1, data[0][1], M1, 1)
plot1(x, M2, ve2, f2, data[0][2], ve2, 2)
plot1(x, M3, ve3, f3, data[0][3], f3, 3)
Na wykresie 1.1 widzimy zmianę drogi przebytej przez rakietę dla różnych jej mas. Im mniejsza masa rakiety, tym droga przebyta przez nią jest większa.
Wykres 2.1 pokazuje jak zmienia się droga dla różnych wartości prędkości spalin. Im większe są te prędkości, tym większą drogę pokonuje rakieta w danym czasie.
Wykres 3.1 przedstawia zmianę drogi dla różnych wartości tempa spalania paliwa. Dla większego tempa spalania, droga przebyta w danym czasie również ma większą wartość.
Droga przebyta przez rakietę w czasie rośnie wykładniczo.
Wykres 1.2 przedstawia jak zmienia się prędkość rakiety dla różnych jej mas. Im mniejsza masa, tym prędkość rakiety wzrasta szybciej w danym czasie.
Na wykresie 2.2 widzimy zmianę prędkości rakiety w zależności od prędkości spalin rakiety. Im prędkość spalin jest większa, tym większą prędkość osiąga rakieta w danym czasie.
Wykres 3.2 pokazuje prędkości rakiety dla różnych wartości tempa spalania paliwa. Jeżeli tempo to zwiększa się, prędkość rakiety również się zwiększa.
Prędkość rakiety w czasie zmienia się logarytmicznie.
Wykres 1.3 przedstawia wartość przyspieszenia rakiety dla różnych jej mas. Im większa masa rakiety, tym przyspiesznie w czasie jest mniejsze.
Na wykresie 2.3 widzimy wartości przyspiesznia w danym czasie dla różnych prędkości spalin rakiety. Im prędkość spalin jest większa, tym większe jest przyspiesznie rakiety w danym czasie.
Wykres 3.3 pokazuje, że im większe jest tempo spalania paliwa, tym przyspieszenie rakiety jest większe.
Przyspiesznie rakiety w czasie zmienia się hiperbolicznie.
from IPython.display import HTML, display
def display_table(data):
"""Create html table from data table given as argument."""
html = "<table>"
for row in data:
html += "<tr>"
for field in row:
html += "<td><h4>%s</h4><td>"%(field)
html += "</tr>"
html += "</table>"
display(HTML(html))
if __name__ == "__main__":
M = [572878, 231420, 28178, 5178]
ve = [2700, 3200, 4400, 2800]
f = [432, 1136, 23.1, 52]
t = [95, 250, 900, 85]
engines = ["SRB", "Atlas V", "Centaur", "Star 48B"]
data = [["silnik", "masa rakiety [kg]", "prędkość spalin [m/s]", "tempo spalania [kg/s]", "czas działania stopnia [s]"]]
for i in range(4):
data.append([engines[i], M[i], ve[i], f[i], t[i]])
display_table(data)
import math
import matplotlib.pyplot as plt
import numpy as np
def velocity_of_multistage_rocket(i, M, ve, f, t):
"""Count multistage rocket velocity.
Args:
i [int]: time from launching rocket (in seconds),
M [list]: masses of rocket with propellant in every stage,
ve [list]: exhaust velocities in every stage,
f [list]: masses of ejected gas per second in every stage,
t [list]: times of working each stage.
Return:
[int]: multistage rocket velocity in i.
"""
if i <= t[0]:
v = (ve[0]+ve[1])/2*math.log(M[0]/(M[0]-5*f[0]*i-f[1]*i))
return v
elif i <= t[1]:
v = ve[1]*math.log(M[1]/(M[1]-f[1]*(i-t[0]))) + (ve[0]+ve[1])/2*math.log(M[0]/(M[0]-5*f[0]*t[0]-f[1]*t[0]))
if i == t[1]:
print("Prędkość rakiety wielostopniowej po pierwszym etapie to {} km/s.".format(round(v/1000,2)))
return v
elif i <= sum(t[1:3]):
v = ve[2]*math.log(M[2]/(M[2]-f[2]*(i-t[1]))) + ve[1]*math.log(M[1]/(M[1]-f[1]*(t[1]-t[0]))) + (ve[0]+ve[1])/2*math.log(M[0]/(M[0]-5*f[0]*t[0]-f[1]*t[0]))
if i == sum(t[1:3]):
print("Prędkość rakiety wielostopniowej po drugim etapie to {} km/s.".format(round(v/1000,2)))
return v
elif i <= sum(t[1:4]):
v = ve[3]*math.log(M[3]/(M[3]-f[3]*(i-sum(t[1:3])))) + ve[2]*math.log(M[2]/(M[2]-f[2]*t[2])) + ve[1]*math.log(M[1]/(M[1]-f[1]*(t[1]-t[0]))) + (ve[0]+ve[1])/2*math.log(M[0]/(M[0]-5*f[0]*t[0]-f[1]*t[0]))
if i == sum(t[1:4]):
print("Ostateczna prędkość rakiety wielostopniowej to {} km/s.".format(round(v/1000,2)))
return v
def plot2(x, M, ve, f, t, M2, ve2, f2):
"""Draw velocity vs. time plot, acceleration vs. time plot and distance vs. time plot for multistage rocket and onestage rocket.
Args:
x [numpy.linspace]: time interval (in seconds),
M [list]: masses of multistage rocket with propellant in every stage,
ve [list]: exhaust velocities in every stage multistage rocket,
f [list]: masses of ejected gas per second in every stage multistage rocket,
t [list]: times of working each stage multistage rocket,
M2 [int]: mass of onestage rocket with propellant,
ve2 [int]: exhaust velocity onestage rocket,
f2 [int]: mass of ejected gas per second onestage rocket.
"""
fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, sharex=True, figsize=(19, 5.5))
multistage_rocket_velocities = [velocity_of_multistage_rocket(i, M, ve, f, t) for i in x]
onestage_rocket_velocities = [velocity(i, M2, ve2, f2) for i in x]
ax0.set_title('4.1.Zależność drogi pokonanej przez rakietę od czasu lotu', fontsize=13)
for list_of_velocities in [onestage_rocket_velocities, multistage_rocket_velocities]:
ax0.plot(x, distance(list_of_velocities))
ax0.set_xlabel('czas [s]', fontsize=10)
ax0.set_ylabel('droga [km]', fontsize=10)
ax0.legend(["Rakieta jednostopniowa", "Rakieta wielostopniowa"])
ax1.set_title('4.2.Zależność prędkości rakiety od czasu lotu', fontsize=13)
for list_of_velocities in [onestage_rocket_velocities, multistage_rocket_velocities]:
ax1.plot(x, list_of_velocities)
ax1.set_xlabel('czas [s]', fontsize=10)
ax1.set_ylabel('prędkość [m/s]', fontsize=10)
ax1.legend(["Rakieta jednostopniowa", "Rakieta wielostopniowa"])
ax2.set_title('4.3.Zależność przyspieszenia rakiety od czasu lotu', fontsize=13)
ax2.plot(x[1:len(x)], acceleration(onestage_rocket_velocities))
ax2.scatter(x[1:len(x)], acceleration(multistage_rocket_velocities), marker='.', color='orange')
ax2.set_xlabel('czas [s]', fontsize=10)
ax2.set_ylabel('przyspieszenie [m/s²]', fontsize=10)
ax2.legend(["Rakieta jednostopniowa", "Rakieta wielostopniowa"])
fig.suptitle('Wykresy zależności - porównanie rakiety wielostopniowej i jednostopniowej', fontsize=16)
plt.show()
if __name__ == "__main__":
x = np.linspace(0, sum(t[1:len(t)]), sum(t[1:len(t)])+1)
plot2(x, M, ve, f, t, 2916500, 34000, 1000)
W rakiecie wielostopniowej w zależności od tego, jaki silnik działa w danym momencie, prędkość rośnie wolniej lub szybciej. Przyspieszenie jest pochodną prędkości, więc jeśli tempo wzrostu prędkości zmniejsza się, przyspieszenie rakiety również maleje.
Końcowa prędkość rakiety oraz droga pokonana przez nią dla rakiety wielostopniowej i jednostopniowej mają zbliżone wartości, jednak warto zauważyć, że masa rakiety wielostopniowej w tym przypadku jest ponad 5 razy mniejsza. Wynika z tego, że rakieta wielostopniowa zostawia mniej śmieci w kosmosie, co jest bardziej ekologicznym rozwiązaniem. ☀️🌱
Równanie opisujące siłę wypadkową działającą na rakietę, przy uwzględnienieniu, że: $$F = F_c - mg - kv^2,$$ gdzie:
Powyższe równanie możemy również zapisać w następującej formie: $$m\frac{dv}{dt} + \frac{dm}{dt}v = F_{c} - mg - kv^2.$$
Niech masa rakiety będzie równa: $$m = m_0 - \gamma t,$$ gdzie:
Podstawiamy $m$ oraz $\frac{dm}{dt}$ do powyższego równania i otrzymujemy: $$(m_0 - \gamma t)\frac{dv}{dt} -\gamma v = F_c - (m_{0} -\gamma t)g - kv^2.$$
Aby znależć numeryczne rozwiązanie równania: $(m_0 - \gamma t)\frac{dv}{dt} -\gamma v = F_c - (m_{0} -\gamma t)g - kv^2$, korzystamy z metody Eulera.
Najpierw przekształcamy równanie do formy: $$\frac{dv}{dt} = \frac{F_c - (m_{0} -\gamma t)g - kv^2 + \gamma v}{m_0 - \gamma t}.$$
Stosujemy metodę Eulera do naszego równania.
Pochodną przybliżamy różnicą skończoną: $$ v'(t_n) \approx \frac{v(t_n + \Delta t)- v(t_n)}{\Delta t} ,$$ gdzie:
przy warunku początkowym: $v(0) = v_0$.
Niech $t_n + \Delta t = t_{n+ 1}$.
Jeżeli $v_n$ oznacza numeryczne przybliżenie do $v(t_n)$ otrzymujemy: $$v_{n + 1} = v_n + \Delta t \frac{dv}{dt}.$$
Podstawimay $\frac{dv}{dt}$ i otrzymujemy: $$v_{n + 1} = v_n + \frac{F_c - (m_{0} -\gamma t_n)g - kv_n^2 + \gamma v_n}{m_0 - \gamma t_n}\Delta t.$$
Użyjemy wzoru, aby przybliżyć numeryczne rozwiązanie równania dla konkretnej rakiety.
Właściwości rakiety:
$F_c = 34000000$ N;
$m_0 = 2500000$ kg;
$\gamma = 1000 \frac{kg}{s}$;
$g = 9,8 \frac{m}{s^2}$;
$k = 0.4 \frac{kg}{m}$;
Niech $\Delta t = 1 s$.
Warunek początkowy: $v(0) = 0$.
Podstawimy dane do równania: $$v_{n + 1} = v_n + \frac{34000000 - (2500000 - 1000t_n)9,8 - 0.4v_n^2 + 1000v_n}{2500000 - 1000t_n} = \frac{9500000 + 98000t_n - 0.4v_n^2 + 1000v_n}{2500000 - 1000t_n }.$$
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [19, 5.5]
from scipy.integrate import odeint
import numpy as np
def euler_velocity(v_prev, max_time):
"""Count rocket velocity using Euler's methode.
Args:
v_prev [int]: initial rocket velocity,
max_time [int]: maximum time value.
Return:
[list]: list of rocket velocities changing every second.
"""
velocities = [v_prev]
for t in range(1, max_time):
v_next = v_prev + ((9500000) + 98000*t - 0.4*(v_prev**2) + 1000*v_prev)/(2500000 - 1000*t)
velocities.append(v_next)
v_prev = v_next
return velocities
def diff_velocity(v, t):
"""Count diffrential value of v at t time."""
dvdt = (9500000 + 98000*t - 0.4*(v**2) + 1000*t)/(2500000 - 1000*t)
return dvdt
def plot3(time_range, solve, euler):
"""Draw velocity vs. time plot.
Args:
time_range [numpy.linspace]: time interval (in seconds),
solve [list]: list of velocities for diffrential equation,
euler [list]: list of velocities for Euler's methode.
"""
plt.plot(time_range, solve)
plt.plot(time_range, euler)
plt.title("5.Wykres zależności prędkości rakiety od czasu lotu", fontsize=16)
plt.xlabel("czas [s]", fontsize=12)
plt.ylabel("prędkość [m/s]", fontsize=12)
plt.legend(["Funkcja Odeint wbudowana w Pythonie", "Metoda Eulera"], fontsize=12)
plt.show()
if __name__=="__main__":
v0 = 0.0
max_time = 600
time_range = np.linspace(0, 600, 600)
solve = odeint(diff_velocity, v0, time_range)
euler = euler_velocity(v0, max_time)
plot3(time_range, solve, euler)
Wykres pokazuje, że metoda Eulera daje dokładne przybliżenia prędkości rakiety dla czasu lotu od 0 do 3 minut. W późniejszym przedziale czasu przybliżanie prędkości tą metodą jest coraz mniej dokładne.
import math
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [19, 5.5]
import numpy as np
def vel(v_init, max_time):
"""Count rocket velocity.
Args:
v_init [int]: initial rocket velocity,
max_time [int]: maximum time value.
Return:
[list]: list of rocket velocities changing every second.
"""
velocities = [v_init]
for t in range(1, max_time):
v = 34000 * math.log(2500000/(2500000 - 1000*t))
velocities.append(v)
return velocities
def plot4(time_range, velocities_tab, max_time):
"""Draw velocity vs. time plot.
Args:
velocity_tab [list]: list of rocket velocities,
time_range [numpy.linspace]: time interval (in seconds),
max_time [int]: maximum time value.
"""
plt.plot(time_range, velocities_tab[0])
plt.plot(time_range, velocities_tab[1])
plt.title("6.Zależność prędkości rakiety od czasu lotu", fontsize=16)
plt.legend(['wartość prędkości, gdy nie działają zewnętrzne siły', 'wartośc prędkości, gdy działa siła grawitacji i oporu atmosferycznego'], fontsize=12)
plt.xlim(0, max_time)
plt.xlabel('czas [s]', fontsize=12)
plt.ylabel('prędkość rakiety [m/s]', fontsize=12)
plt.show()
if __name__=="__main__":
time_range = np.linspace(0, 180, 180)
v0 = 0
max_time = 180
plot4(time_range, [vel(v0, max_time), euler_velocity(v0, max_time)], max_time)
Wykres pokazuje prędkość rakiety w przypadku, kiedy nie działałyby na nią żadne zewnętrzne siły oraz w przypadku, gdy działa na nią siła grawitacji i siła oporu atmosferycznego.
Przy działaniu sił zewnętrznych prędkość rakiety w danym czasie jest o wiele mniejsza, niż kiedy nie działają na nią żadne siły.
W naszym projekcie wykorzystujemy równanie Ciołkowskiego, które jest bazą przy modelowaniu ruchu rakiet oraz pojazdów zachowujących się jak rakiety. Jest to podstawowy wzór stosowany w tych modelach. Oczywiście, projektując rakietę, uwzględnia się również inne czynniki, takie jak działanie sił oporu, przyciąganie grawitacyjne, ciśnienie atmosferyczne czy parametry silnika rakiety. W naszej pracy przedstawiamy również model rakiety z uwględnieniem siły grawitacji i oporu atmosferycznego oraz porównujemy jej zachowanie w przypadku, kiedy nie działają na nią żadne zewnętrzne siły.
Przedstawiamy również model rakiety wielostopniowej, którą porównujemy z rakietą jednostopniową. Ponadto sprawdzamy dokładność metody numerycznej Eulera z uwzględnieniem działania na rakietę zewnętrznych sił, co prezentujemy na wykresie.
Wynika z tego, że model lotu rakiety może być bardzo złożony i w rzeczywistości konstruktorzy rakiet muszą uwzględnić wiele dodatkowych czynników, które wpływają na jej lot.