LOT RAKIETY 🚀

„Jeśli potrafisz wynieść swój statek na orbitę, jesteś w połowie drogi dokądkolwiek.”

Robert Heinlein

Autorzy:

Aleksandra Janczewska,
Katarzyna Turbańska

Wstęp

Od prochu strzelniczego do misji kosmicznych

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.

Jak właściwie działa rakieta?

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.

Definicje i wzory

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]}$:

  • możliwość skorzystania z innego typu silnika - lepiej dopasowanego do warunków jego pracy - dla każdego stopnia. W „dolnych” członach stosuje się silniki dostosowane do ciśnienia atmosferycznego, w „górnych” natomiast - silniki zaprojektowane do działania w warunkach bliskich próżni;
  • pozbywanie się zbędnej masy w trakcie lotu;
  • według Konstantina Ciołkowskiego rakieta jednostopniowa, ze względu na zbyt duży ciężar, nie jest w stanie osiągnąć orbity ziemskiej.

Wady rakiet wielostopniowych$^{[7]}$:

  • zwiększenie stopnia skomplikowania konstrukcji, co skutkuje większym zagrożenie awariami - nieprawidłowe oddzielenie któregoś ze stopni, kolizje poszczególnych członów bądź niepoprawny zapłon kolejnego silnika;
  • konieczność przenoszenia dodatkowego ciężaru w postaci silników niewykorzystywanych w większości faz lotu.

Przyspieszenie rakiety z uwględnieniem siły grawitacji$^{[5]}$$^{[8]}$

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:

  • im większa jest wartość prędkości wyrzucanego paliwa względem rakiety ($v_{e}$), tym większe jest przyspieszenie rakiety. Granica $v_{e}$ dla konwencjonalnych układów napędowych na gorący gaz wynosi około $2,5 \cdot 10^3\frac{m}{s}$;
  • iloczyn $(\frac{dm}{dt})v_{e}$ to siła ciągu rakiety. Im szybciej rakieta spala paliwo, tym większa jest ta siła, z czego wynika, że przyspieszenie rakiety jest większe;
  • masa rakiety $m$ drastycznie spada podczas lotu, ponieważ jej większość to paliwo. Przyspieszenie rośnie w sposób ciągły, osiągając maksimum tuż przed wyczerpaniem się paliwa.

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 Ciołkowskiego$^{[5]}$

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:

  • $\Delta v$ - końcowa (maksymalna) prędkość pojazdu;
  • $m_{0}$ - całkowita masa początkowa pojazdu wraz z paliwem;
  • $m_{1}$ - masa końcowa pojazdu bez paliwa;
  • $v_{e}$ - prędkość wydzielanego czynnika roboczego (spalonego paliwa) względem pojazdu.

Przykład$^{[8]}$:

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.

Wyprowadzenie wzoru Ciołkowskiego$^{[9]}$

Rysunek poglądowy

Wiemy, że prędkość spalin względem rakiety ($v_{e}$) dana jest zależnością: $$v_{e} = v - v_{s},$$ gdzie:

  • $v_{s}$ - prędkość spalin opuszczających silnik względem obserwatora (na Ziemi);
  • $v$ - prędkość chwilowa rakiety względem obserwatora (na Ziemi).

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:

  • $\Delta v$ - prędkość końcowa rakiety;
  • $m_{0}$ - masa początkowa rakiety z paliwem;
  • $m_{1}$ - masa końcowa rakiety bez paliwa.

Ostatecznie otrzymujemy: \begin{equation} \Delta v = v_{e}\ln{\frac{m_{0}}{m_{1}}}. \end{equation}

Wykresy zależności bez uwzględniania sił zewnętrznych działających na rakietę

Wzór na zmianę prędkości rakiety w czasie: $$v = v_e\ln\left( \frac{m}{m - \frac{dm}{dt}t} \right),$$ gdzie:

  • $t$ - czas od początku startu rakiety;
  • $m$ - masa początkowa rakiety z paliwem.

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.

Tabele właściwości rakiet jednostopniowych (dane do wykresów)

In [43]:
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("---------------------------------------------------------------------")

TABELA 1

masa rakiety [kg]

prędkość spalin [m/s]

tempo spalania [kg/s]

rakieta 1

2500000

34000

1000

rakieta 2

2916500

34000

1000

rakieta 3

3100000

34000

1000

rakieta 4

3256500

34000

1000

---------------------------------------------------------------------

TABELA 2

masa rakiety [kg]

prędkość spalin [m/s]

tempo spalania [kg/s]

rakieta 1

2916500

31000

1000

rakieta 2

2916500

32000

1000

rakieta 3

2916500

33000

1000

rakieta 4

2916500

34000

1000

---------------------------------------------------------------------

TABELA 3

masa rakiety [kg]

prędkość spalin [m/s]

tempo spalania [kg/s]

rakieta 1

2916500

34000

1000

rakieta 2

2916500

34000

1200

rakieta 3

2916500

34000

1400

rakieta 4

2916500

34000

1600

---------------------------------------------------------------------
In [50]:
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)

Opis wykresów

1. Droga, którą pokonuje rakieta w danym czasie.

  • 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.

2. Prędkość, którą osiąga rakieta w danym czasie.

  • 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.

3. Przyspieszenie, które osiąga rakieta w danym czasie.

  • 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.

Wielostopniowa rakieta nośna sondy New Horizons$^{[10]}$

Silniki SRB oraz Atlas V zaczynają pracę jednocześnie. Atlas V pracuje dalej po odrzuceniu SRB.

Tabela właściwości wielostopniowej rakiety nośnej sondy New Horizons (dane do wykresów)

In [51]:
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)

silnik

masa rakiety [kg]

prędkość spalin [m/s]

tempo spalania [kg/s]

czas działania stopnia [s]

SRB

572878

2700

432

95

Atlas V

231420

3200

1136

250

Centaur

28178

4400

23.1

900

Star 48B

5178

2800

52

85

In [52]:
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)
Prędkość rakiety wielostopniowej po pierwszym etapie to 6.91 km/s.
Prędkość rakiety wielostopniowej po drugim etapie to 12.8 km/s.
Ostateczna prędkość rakiety wielostopniowej to 18.18 km/s.

Opis wykresów

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 rakiety jednostopniowej z uwzględnieniem siły grawitacji i oporu atmosferycznego

Równanie opisujące siłę wypadkową działającą na rakietę, przy uwzględnienieniu, że: $$F = F_c - mg - kv^2,$$ gdzie:

  • $F = \frac{d(mv)}{dt}$ - siła wypadkowa;
  • $F_c = const$ - siła ciągu rakiety;
  • $g = const$ - przyspieszenie grawitacyjne;
  • $m$ - masa rakiety zmieniająca się w czasie;
  • $k$ - współczynnik oporu atmosferycznego;
  • $v$ - prędkość rakiety zmieniająca się w czasie.

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:

  • $m_0$ - masa początkowa rakiety (masa rakiety z paliwem);
  • $\gamma = \frac{dm}{dt}$ - zmiana masy rakiety w czasie, to znaczy wartość tempa spalania paliwa;
  • $t$ - czas lotu rakiety.

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.$$

Metoda Eulera

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:

  • $\Delta t$ - przyrost czasu,

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 }.$$

Wykres pokazujący dokładność metody numerycznej Eulera

In [63]:
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)
    

Analiza wykresu

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.

Wykres zależności prędkości od czasu w zależności od sił działających na rakietę

In [48]:
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)

Analiza wykresu

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.

Podsumowanie

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.