Zadanie

Pamätáte sa na piráta Francisa, ako so svojou loďou plával po plochej Zemi? Nepamätáte? No čo už. Ako si tak plával cez Pacifik, odchytili ho aj s loďou mimozemšťania pomocou ťažného lúča a uskladnili za účelom neskoršej analýzy. Francisova loď sa teraz voľne vznáša v bezváhovom stave vnútri gigantickej kozmickej lode. Kam len Francis dovidí, rozprestiera sa len obyčajný vzduch s teplotou \(\SI{0}{\celsius}\) pri tlaku \(\SI{100}{\kilo\pascal}\).

V rámci protestu proti tomuto uväzneniu Francis prikáže svojim námorníkom, aby z dela vypálili olovenú guľu s priemerom \(\SI{10}{\centi\metre}\). Guľa vyletí z hlavne rýchlosťou \(\SI{100}{\metre\per\second}\). Ako ďaleko doletí, než ju odpor vzduchu zastaví a ako dlho to potrvá, ak odporová sila závisí od rýchlosti ako \(F = -\frac{1}{2} C_d S \rho_\mathrm{vzduch} v^2\)?

Najprv uvažujte konštantný odporový koeficient pre guľu, \(C_d = \num{0.47}\). Námorní delostrelci však dobre vedia, že to také jednoduché nie je. V šestnástom storočí možno odporový koeficient bol konštantný, doba však pokročila a dnes už závisí aj od Reynoldsovho čísla. Konkrétne takto, podľa Morrison, 2016: \[ C_d(\mathrm{Re}) = \frac{24}{\mathrm{Re}} + \frac{ \num{2.6} \left(\frac{\mathrm{Re}}{5}\right) }{ 1 + \left(\frac{\mathrm{Re}}{5}\right)^{\num{1.52}} } + \frac{ \num{0.411} \left(\frac{\mathrm{Re}}{263000}\right)^{\num{-7.94}} }{ 1 + \left(\frac{\mathrm{Re}}{263000}\right)^{-8} } + \frac{ \num{0.25} \left(\frac{\mathrm{Re}}{\num{e6}}\right) }{ 1 + \left(\frac{\mathrm{Re}}{\num{e6}}\right) }. \]

Uvažujte, že po vypálení gule zostane Francisova loď v pokoji. Pri riešení úlohy odporúčame programovať alebo použiť nejaký tabuľkový procesor.

Úloha vyzerá celkom priamočiaro – a aj je, akurát si treba dať pozor a nenechať sa oklamať čiastkovým úspechom. Simulovať môžeme buď v nejakom programovacom jazyku, alebo sa uspokojíme s nejakým tabuľkovým procesorom, napríklad Excelom. Simulácia bude v oboch prípadoch rovnaká, akurát za \(C_d\) prvý raz dosadíme konštantu, druhý raz si doň musíme v každom kroku priradiť výsledok pomerne zložitého výpočtu. Ten však počítač urobí za nás, stačí mu raz povedať, ako sa to robí.

Keďže kroky si môžeme voliť ako chceme1, postačí nám jednoduchá Eulerova metóda: v každom kroku si analyticky spočítame zrýchlenie \(a\). Tu našťastie závisí iba od rýchlosti ako \[ a(v) = \frac{F(v)}{m} = -\frac{C_d S \rho_\mathrm{a} v^2}{2 m} = -\frac{C_d \pi r^2 p M_{\mathrm{a}} v^2}{2 m R T}, \qquad(1)\] kde sme hustotu vyjadrili zo stavovej rovnice ideálneho plynu, pričom sa však nič hrozné nestane, ak za ňu dosadíme tabuľkovú hodnotu.

Teraz v každom kroku rýchlosť zväčšíme o \(a \Delta t\)2, polohu zväčšíme o \(v \Delta t\) a tak až donekonečna – alebo teda kým neuvidíme, že sa delová guľa zastavila, resp. jej rýchlosť klesla pod nejakú kritickú hodnotu, povedzme \(\SI{1}{\micro\metre\per\second}\). Teda máme predpis \[ \begin{aligned} x_0 &= \SI{0}{\metre}, \qquad &v_0 &= \SI{100}{\metre\per\second}, \\ x_{n + 1} &\leftarrow x_n + v_n \cdot \Delta t, \qquad &v_{n + 1} &\leftarrow v_n + a(v_n) \cdot \Delta t \end{aligned} \qquad(2)\] a s tým si už program alebo tabuľkový procesor hravo poradia.

Jednoduché? Vyzerá to tak. Lenže pri numerickom riešení diferenciálnych rovníc často natrafíme na nepríjemný jav. Presnosť totiž závisí od veľkosti kroku.

  • Ak zvolíme priveľký krok, pri hodnotách, kde sa zrýchlenie prudko mení (v našom prípade na začiatku pohybu, keď je odporová sila veľká), stratíme presnosť. Ak by sme zvolili veľmi veľký krok, guľa by v našom výpočte zastavila prakticky okamžite, alebo sa dokonca začala pohybovať naspäť.
  • Ak zvolíme primalý krok, výpočet síce bude veľmi presný… akurát jeho konca sa nemusíme dožiť.

Na tejto rovnici je trochu nepríjemné práve to, že na dosiahnutie dostatočnej presnosti potrebujeme veľmi malý krok3. Ak sa uspokojíme s prvým tipom, ktorý zvolíme, určite nedostaneme dobrý výsledok. Ak \(\Delta t\) trochu zväčšíme, zrazu nám vyjde niečo úplne iné. Aby nám malá zmena \(\Delta t\) dala iba malú zmenu výsledku, musíme použiť krok natoľko malý, že nás pred ukončením výpočtu zastihne dôchodkový vek. Rozumná presnosť sa dá dosiahnuť pre \(\Delta t = \SI{0.001}{\second}\). Pri sekundovom kroku už tratíme veľa presnosti na začiatku pohybu.

S trochou zamyslenia sa ale vieme vynájsť: čo keby sme použili malé časové kroky vtedy, keď sa zrýchlenie prudko mení (teda zo začiatku), a veľké kroky neskôr, keď budú rýchlosť a odporová sila menšie? Napríklad po každom kroku zväčšíme \(\Delta t\) o , čiže ho vynásobíme \(\num{1.01}\). Takto vieme dosiahnuť obstojnú presnosť a zároveň stihnúť termín odovzdania. Rovnica 2 samozrejme prejde na \[ \begin{aligned} x_0 &=& \SI{0}{\metre}, \qquad v_0 &=& \SI{100}{\metre\per\second}, \qquad \Delta t_0 &=& \SI{0.001}{\second}, \\ x_{n + 1} &\leftarrow& x_n + v_n \cdot \Delta t_n, \qquad v_{n + 1} &\leftarrow& v_n + a(v_n) \cdot \Delta t_n, \quad \Delta t_{n + 1} &\leftarrow& \Delta t_n \cdot \num{1.01}. \end{aligned} \qquad(3)\]

Prvý prípad

V prvom prípade by sme mali zistiť, že nech robíme, čo chceme, vzdialenosť nám vždy závisí od počtu krokov, ktoré necháme počítač odsimulovať. Dokonca pekne exponenciálne: ak celkovú dobu simulácie zväčšíme \(e\)-krát, guľa zaletí o skoro presne \(\SI{2500}{\metre}\) ďalej. Guľa sa teda nezastaví nikdy a preletená vzdialenosť môže byť ľubovoľne veľká, ak jej dáme dosť času.

Od vás sme to samozrejme nechceli, ale táto rovnica sa dá vyriešiť aj analyticky: \[ \frac{\mathrm{d} v}{\mathrm{d} t} = -\alpha v^2 \quad\text{s počiatočnou podmienkou}\quad v(0) = \SI{100}{\metre\per\second} \quad\Rightarrow\quad v(t) = \frac{1}{100 \alpha t + 1} \cdot \SI{100}{\metre\per\second}. \]

Pre \(t \to \infty\) síce rýchlosť klesá k nule, ale nie dosť rýchlo. Prejdená dráha je rovná \[ x(t) = \int\limits_{0}^{t} v(\tau) \mathrm{d} \tau = \frac{\ln\left(100 \alpha t + 1\right)}{\alpha}, \] a teda rastie síce veľmi pomaly, ale nad všetky medze.

Druhý prípad

Popísaný vzťah je nechutný, ale zato veľmi dobre sedí s experimentálnymi dátami. Ani v tomto prípade sa guľa nikdy úplne nezastaví, pokles rýchlosti je ale dostatočne rýchly na to, aby celková vzdialenosť bola konečná. Žiadnym nastavením presnosti alebo dĺžky výpočtu nedosiahneme, aby preletela viac ako \(\SI{32.639}{\kilo\metre}\).

Aj keď kompletné analytické riešenie je v tomto prípade nedosiahnuteľné4, pre malé \(\mathrm{Re}\) sa funkcia \(C_d\) rovná zhruba \(\frac{24}{\mathrm{Re}}\) a odporová sila sa približuje k Stokesovmu odporu \[ F = - \underbrace{6 \pi \mu r}_{\beta} v. \]

Ako guľa brzdí, skôr či neskôr musí spomaliť natoľko, aby viskozita vzduchu začala hrať dôležitú rolu. Pohyb je potom popísaný rovnicou \[ \frac{\mathrm{d} v}{\mathrm{d} t} = -\beta v \quad\Rightarrow\quad v(t) = v(t_0) \mathrm{e}^{-\beta t}, \]

Rýchlosť klesá k nule, aj keď ani teraz ju nikdy nedosiahne. Prejdená dráha je však rovná \[ x(\infty) = \int\limits_{t_0}^{\infty} v(\tau) \mathrm{d} \tau = \frac{v(t_0)}{\beta}, \] čo celkom veselo konverguje ku konečnej hodnote.

Komentár k riešeniam

Úloha mala dve podstatné podčasti: správne naprogramovať samotné riešenie, a potom si zvoliť správnu dĺžku kroku. Pravda, nešlo o to, či nájdete konkrétnu hodnotu alebo predpis, ktorý dá ako-tak správny výsledok, ale či si všimnete, čo to robí.


  1. Vychádzame z trúfalého predpokladu, že ste nezačali riešiť o 23:30 v posledný deň pred termínom.↩︎

  2. Zrýchlenie \(a\) bude záporné, takže ju v skutočnosti budeme znižovať.↩︎

  3. Áno, v prvej poznámke pod čiarou som klamal.↩︎

  4. Museli by sme totiž integrovať čosi, čo závisí aj od tej príšernej funkcie zo zadania.↩︎

Diskusia

Tu môžte voľne diskutovať o riešení, deliť sa o svoje kusy kódu a podobne.

Pre pridávanie komentárov sa musíš prihlásiť.