Keďže simulovať pohyb úlomkov je jednoduchšie ako riešiť nehomogénne diferenciálne rovnice, je jasné, ktorý spôsob riešenia si vyberieme. Bude to úplne najjednoduchšia Eulerova metóda. Tá nám hovorí, že ak poznáme v čase $t$ zrýchlenie úlomka $\vec{a}(t)$, tak ak úlomok mal v danom čase rýchlosť $\vec{v}(t)$ a polohu $\vec{r}(t)$, tak o malý čas $\FDiff t$ neskôr už mal rýchlosť $$ \vec{v}(t + \FDiff t) = \vec{v}(t) + \vec{a}(t)\FDiff t $$ a polohu $$ \vec{r}(t + \FDiff t) = \vec{r}(t) + \vec{v}(t)\FDiff t. $$ Na malom časovom úseku sme teda nejaký komplikovaný pohyb nahradili rovnomerným priamočiarym. A potom ak chceme poznať polohu a rýchlosť v čase $2\FDiff t$, tak budeme vychádzať z hodnôt v čase $\FDiff t$. Ak zvolíme časový krok dostatočne malý, tak naše numerické riešenie dosiahne nejakú presnosť, ktorú požadujeme. A to je všetko. Jednoduché, však?

Doteraz sme veličiny zapisovali všeobecne ako vektory, ale teraz si ich rozdelíme na zložky. Počiatočná poloha všetkých úlomkov je rovnaká. Danú máme výšku $h$ a vodorovnú súradnicu zvoľme nulovú, čiže $$ x(0) = 0 \QQText{a} y(0) = h. $$ Výbuchom dostanú všetky úlomky hybnosť $\FDiff p$. Každý má rovnakú hmotnosť $m$, čiže úlomky výbuch urýchlil na rýchlosť $\FDiff p/m$ voči ponorke. Tá išla rýchlosťou $u$ voči vode, ktorá ide rýchlosťou $v$ voči dnu. Preto počiatočné zložky rýchlosti sú $$ v_x(0) = u + v + \frac{\FDiff p}{m}\cos\alpha \QQText{a} v_y(0) = \frac{\FDiff p}{m}\sin\alpha, $$ kde $\alpha$ je uhol, pod ktorým úlomok vyletel, rátaný v kladnom smere od smeru rýchlosti ponorky. Ohľadom zrýchlenia použijeme Stokesov zákon $\vec{F} = -6\pi\mu R \vec{v}$, ktorý nám umožňuje vypočítať odporovú silu pôsobiacu na guľu polomeru $R$ pohybujúcu sa rýchlosťou $v$ voči kvapaline s dynamickou viskozitou $\mu$. Zrýchlenie je $\vec{a} = \vec{F}/m$ a guľa s hmotnosťou $m$ má polomer $R = \sqrt[3]{3m/(4\pi\rho)}$, kde $\rho$ je hustota gule. Okrem odporovej sily pôsobia ešte gravitačná a vztlaková v zvislom smere. Dokopy teda máme vo vodorovnom smere $$ a_x(t) = \frac{-6\pi\mu\sqrt[3]{\frac{3m}{4\pi\rho}}(v_x(t) - v)}{m}. $$ Rozdiel $v_x(t) - v$ je vodorovná zložka relatívnej rýchlosti úlomku voči vode. Pre zvislý smer máme $$ a_y(t) = \frac{-6\pi\mu\sqrt[3]{\frac{3m}{4\pi\rho}} v_y(t)}{m} - g + \frac{\rho_\mathrm{v}}{\rho}g, $$ kde posledný člen je pochádza zo vztlakovej sily.

Časové kroky pre každý úlomok vykonávame, kým platí $y > 0$, teda že je nad dnom. Ak táto podmienka prestane platiť, tak pre aktuálne hodnoty vypočítame vzdialenosť od miesta výbuchu Pytagorovou vetou $d = \sqrt{x^2 + h^2}$.

Pre hodnoty $h = \qty{3}{\metre}$, $v = \qty{6}{\metre\per\second}$, $u = \qty{2}{\metre\per\second}$, $\mu = \qty{8.9e-4}{\pascal\second}$ (voda), $m = \qty{1}{\kilo\gram}$, $\rho = \qty{7850}{\kilo\gram\per\metre\cubed}$ (oceľ), $\rho_\mathrm{v} = \qty{1000}{\kilo\gram\per\metre\cubed}$, $\FDiff t = \qty{1e-4}{\second}$, $g = \qty{10}{\metre\per\second\squared}$ a $\FDiff p = \qty{10}{\kilo\gram\metre\per\second}$ dostávame, že ak si celých $\ang{360}$ rovnomerne rozdelíme na 1000 uhlov $\alpha$, tak najďalej doletel úlomok, ktorý bol vyletel pod uhlom $\ang{52.61}$, a to konkrétne do vzdialenosti $\qty{30.28}{\metre}$.

Na obrázku 1 možno videť dráhy ôsmich vybraných úlomkov s rôznymi uhlami $\alpha$. Aby bolo vidno unášavé pôsobenie vody, použili sme preň tisícnásobok dynamickej viskozity vody.

obrázok 1: Dráhy ôsmich vybraných úlomokov

Kód riešiaci túto úlohu takisto prikladáme.

import numpy as np                 # kniznica s matematickymi funkciami

h = 3                              # vyska nad dnom
v = 6                              # rychlost vody voci dnu
u = 2                              # rychlost ponorky voci vode
visk = 8.9e-4                      # dynamicka viskozita vody
m = 1                              # hmotnost guliciek
rho = 7850                         # hustota guliciek
R = np.cbrt(3*m / (4*np.pi*rho))   # polomer guliciek
rhov = 1000                        # hustota vody
dt = 0.0001                        # casovy krok
g  = 10                            # gravitacne zrychlenie
dp = 10                            # udelena hybnost
vp = dp/m                          # udelena rychlost

vzdialenosti = []                  # zoznam na ukladanie vzdialenosti od miesta vybuchu
uhly = np.linspace(0, 360, 1000)   # pocitame pre 1000 uhlov medzi 0 a 360 stupnami
for a in uhly:
    # pociatocna poloha
    x = 0
    y = h
    # pociatocna rychlost
    vx = u + v + vp*np.cos(np.radians(a))
    vy = 0     + vp*np.sin(np.radians(a))

    # simulujeme, dokym je gulicka nad dnom
    while(y > 0):
        ax = -6*np.pi*visk*R*(vx - v)/m
        ay = -6*np.pi*visk*R*vy/m - g + rhov/rho*g
        vx = vx + ax*dt
        vy = vy + ay*dt
        x  =  x + vx*dt
        y  =  y + vy*dt

    vzdialenosti.append(np.sqrt(x**2 + h**2)) # ulozime vzdialenost od miesta vybuchu

index = np.argmax(vzdialenosti)    # index najvacsej hodnoty v zozname vzdialenosti
print(f"Najvacsiu vzdialenost {vzdialenosti[index]} m presla gulicka vyletevsia pod uhlom {uhly[index]}°.")