8.2 Drag og 2D-baner numerisk

Bidraget fra luftmotstanden til akselerasjonslinja, for den kvadratiske modellen:

a = ... - (D/m)*abs(v[i])*v[i]

Uttrykket kommer rett fra kraftmodellen $F_D = -D\lvert v\rvert v$ delt på massen. Dragkoeffisienten $D$ er ikke masseuavhengig, så den må deles på $m$ — i motsetning til tyngdeleddet, der massen forkortes.

For den lineære modellen blir leddet - (kv/m)*v[i]. Der trengs ingen absoluttverdi: farten står allerede i første potens, og fortegnet følger av seg selv.

Fordi $v^2$ er positiv uansett fortegn på farten, mens luftmotstanden må skifte fortegn samtidig med bevegelsen.

Sett inn en negativ fart, altså et legeme på vei nedover med $y$ positivt oppover:

Å føre Euler–Cromer-løkka én gang for hver romlig retning, med egne tabeller for hver komponent:

x, y, vx, vy

Alle fire oppdateres i samme gjennomløp av løkka, med samme $dt$ og samme indeks $i$. Rekkefølgen internt er den samme som i én dimensjon: begge fartskomponentene først, deretter begge posisjonskomponentene med de oppdaterte fartene.

Dette er poengfellen sensor er hardest på. En kode som bare tar med én romlig komponent i et todimensjonalt problem, gir maks 2 av 5 poeng, uansett hvor korrekt den ellers er.

$$v = \sqrt{v_x^2 + v_y^2} \qquad \longrightarrow \qquad \texttt{v = np.sqrt(vx[i]**2 + vy[i]**2)}$$

Denne størrelsen er en skalar — den har ingen retning og skrives derfor uten fet markering. Den lagres ikke i noen tabell, siden den regnes på nytt hver gang og ikke trengs etterpå.

Det er nettopp her de to retningene kobles: $v$ inneholder begge komponentene, og $v$ inngår i begge akselerasjonslinjene.

$$F_{D,x} = -D\,v\,v_x, \qquad F_{D,y} = -D\,v\,v_y$$

der $v = \sqrt{v_x^2+v_y^2}$ er tallverdien av farten.

Merk hvor absoluttverdien er blitt av. I én dimensjon skrev vi $-D\lvert v\rvert v$. I to dimensjoner er $v$ per definisjon positiv (det er en lengde), og retningen kommer fra fortegnet på komponenten $v_x$ eller $v_y$. Produktet $v\,v_x$ gjør altså nøyaktig samme jobb som $\lvert v\rvert v$ gjorde i én dimensjon — automatisk riktig fortegn i begge retninger, hele veien.

I kode:

ax = -(D/m)*v*vx[i]
ay = -g - (D/m)*v*vy[i]

$$v_{x,0} = v_0\cos\theta, \qquad v_{y,0} = v_0\sin\theta$$

vx[0] = v0*np.cos(theta)
vy[0] = v0*np.sin(theta)

Vinkelen må være i radianer. np.cos og np.sin regner ikke i grader, og np.radians(38.0) er den trygge omskrivingen. En vinkel oppgitt i grader og satt rett inn er en av de feilene som gir helt plausible, men helt gale tall.

Kontrollen er billig: $v_{x,0}^2 + v_{y,0}^2$ skal gi $v_0^2$ tilbake.

Testen som avslutter løkka når legemet treffer bakken:

if y[i+1] < 0:
    break

Løkka oppdager først at $y$ er blitt negativ, altså i steget etter nedslaget. Med $dt = 0{,}0010\ \text{s}$ og en vannrett fart på $20\ \text{m/s}$ er overskytingen under to centimeter — normalt godt innenfor det du oppgir med tre gjeldende siffer.

Trenger du nedslagspunktet mer nøyaktig, interpolerer du lineært mellom de to siste punktene: brøkdelen $f = y_{n-1}/(y_{n-1}-y_n)$ av det siste steget, og nedslaget ved $x_{n-1} + f\,(x_n - x_{n-1})$.

$$a_x = -\frac{GM\,x}{\left(x^2+y^2\right)^{3/2}}, \qquad a_y = -\frac{GM\,y}{\left(x^2+y^2\right)^{3/2}}$$

I kode:

r = np.sqrt(x[i]**2 + y[i]**2)
ax = -GM*x[i]/r**3
ay = -GM*y[i]/r**3

Legg merke til tredje potens i nevneren. Kraften går som $1/r^2$, og faktoren $x/r$ er komponenten av enhetsvektoren — til sammen $x/r^3$. Å skrive $r^2$ i nevneren og glemme normeringen er en av de vanligste feilene i denne innpakningen.

Strukturen er ellers identisk med kastet: fire tabeller, farten oppdatert først, posisjonen med den oppdaterte farten. Gravitasjon og sentralkraftbevegelse behandles i Del 9, som ikke er skrevet ennå — men koden trenger du ikke vente på, for den er den samme malen.

$$v = \sqrt{v_x^2+v_y^2}$$

Uten drag er $dv_x/dt = 0$ og $dv_y/dt = -g$ to helt separate likninger som kan løses hver for seg. Med drag kan ingen av dem løses uten den andre.

Koblingen er også grunnen til at numerikk fungerer så godt her: en maskin bryr seg ikke om at to likninger er koblet, så lenge begge kan skrives ned.

At oppturen og nedturen ikke er speilbilder av hverandre når luftmotstanden er med.

For ballen i eksempel 2 tar oppturen $1{,}69\ \text{s}$ og nedturen $1{,}86\ \text{s}$, og toppunktet ligger nærmere utkastet enn nedslaget. Nedturen er brattere.

Årsaken er at farten er størst i første del av banen, og dragkraften går som kvadratet av farten. Bremsingen er derfor kraftigst tidlig, og ballen har mistet mye av den vannrette farten før den begynner å falle.

Vinkelen som gir størst rekkevidde, som ikke lenger er $45^\circ$ når luftmotstanden er med.

For ballen i eksempel 2 ligger optimum ved om lag $41$–$42^\circ$. Jo større dragkoeffisient og jo større startfart, desto lavere blir vinkelen.

Maksimumet er dessuten flatt: mellom $38^\circ$ og $45^\circ$ varierer rekkevidden med under én prosent. Det gjør spørsmålet praktisk uinteressant for en kaster, men prinsipielt interessant — det finnes ingen formel, bare en kjøring.

At en ball som kastes og lander i samme høyde, kommer ned med lavere fart enn den ble kastet med, så snart luftmotstand er med.

Uten motstand er de to like, fordi mekanisk energi er bevart og høyden er den samme. Med motstand gjør dragkraften negativt arbeid i begge faser, og differansen tas fra den kinetiske energien.

For ballen i eksempel 2 er startfarten $32{,}0\ \text{m/s}$ og nedslagsfarten $21{,}4\ \text{m/s}$ — et tap på en tredel. Oppgi alltid nedslagsfarten med retning: her skrått nedover, om lag $48^\circ$ under vannrett.

$$x_1 = \frac{v_0^2}{g}\sin 2\theta, \qquad t_1 = \frac{2v_0}{g}\sin\theta, \qquad h = \frac{(v_0\sin\theta)^2}{2g}$$

For eksempel 2 gir koden $101{,}26\ \text{m}$ mot formelens $101{,}28\ \text{m}$ — samsvar på fire gjeldende siffer.

Vær klar over hva kontrollen ikke tester. Med $D = 0$ er hele dragleddet slått av, og en hvilken som helst feil der ville passert uoppdaget. Dragleddet må testes separat, med et endimensjonalt fall der $\tanh$-løsningen er fasit.

$$u_x = v_x - w, \qquad u_y = v_y, \qquad u = \sqrt{u_x^2+u_y^2}$$

Dragleddene blir $-Du\,u_x$ og $-Du\,u_y$, mens tyngden og posisjonsoppdateringen er uendret: legemet flytter seg fortsatt med sin egen fart over bakken.

Vinden virker altså bare gjennom luftmotstanden. I vakuum har den ingen virkning i det hele tatt — og en medvind endrer både flytiden og høyden, ikke bare den vannrette avstanden.

$$F_{D,x} = -k_v v_x, \qquad F_{D,y} = -k_v v_y$$

Her er komponentene ikke koblet: hver likning inneholder bare sin egen fartskomponent, og hver kan løses analytisk for seg. Det er nettopp derfor den lineære modellen er den som brukes når en oppgave vil ha et lukket svar.

Merk kontrasten: i den lineære modellen trengs verken absoluttverdi eller fellesfaktor $v$, fordi farten står i første potens og fortegnet følger av seg selv.

$$-D\,v\,v_x \quad \text{har alltid motsatt fortegn av } v_x$$

Uttrykket gjør derfor nøyaktig samme jobb som $-D\lvert v\rvert v$ gjorde i én dimensjon, og det er den samme mekanismen: størrelsen fra en positiv faktor, retningen fra komponenten.

Skriver du i stedet $-Dv_x^2$, mister du både fortegnsbyttet og koblingen til den andre komponenten — to feil i ett uttrykk.

$$[D] = \text{kg/m}, \qquad [k_v] = \text{kg/s}$$

Koden regner uten å vite noe om enheter, så et bytte gir ingen feilmelding — bare et pent tall som er galt.

Den raskeste kontrollen er terminalfarten: $v_T = \sqrt{mg/D}$ for den kvadratiske modellen og $v_T = mg/k_v$ for den lineære. Får du en terminalfart som er absurd stor eller liten, har du sannsynligvis brukt feil modell.

Strukturen som skiller en 2D-kode fra en 1D-kode:

$$a_x = -\frac{GM\,x}{r^3}, \qquad a_y = -\frac{GM\,y}{r^3}, \qquad r = \sqrt{x^2+y^2}$$

Tredje potens i nevneren, ikke andre. Kraften går som $1/r^2$, og faktoren $x/r$ er $x$-komponenten av enhetsvektoren inn mot sentrum — til sammen $x/r^3$.

Den vanligste feilen er å skrive $r^2$ i nevneren. Symptomet er en bane som er kvalitativt gal på en måte som er lett å se: kraften avtar da altfor sakte med avstanden, og et legeme som skulle unnsluppet, blir fanget.

$$v = \sqrt{\frac{GM}{r_0}}$$

vinkelrett på radien. Avstanden $r$ skal da holde seg konstant gjennom hele simuleringen.

Testen er sterk fordi den treffer tre ting samtidig: begge komponentene (en sirkelbane bruker dem like mye), fortegnene (feil fortegn gir en bane som flyr av gårde), og metoden (med Forward Euler vokser radien jevnt selv i det som skulle vært en perfekt sirkel).

Omløpstiden skal bli $T_p = 2\pi r_0/v$, som gir en fjerde, uavhengig kontroll.

Skillet som avgjør om et problem trenger to tabeller eller fire.

En kloss på et skråplan, en perle på en tråd eller en sykkel i en bakke er endimensjonale: bevegelsen er bundet til én linje, og alle krefter dekomponeres langs og på tvers av den, der tverrkomponentene opphever hverandre.

En ball i lufta, en komet eller en satellitt er todimensjonale: ingenting binder dem til en linje.

Det er antallet frie retninger som avgjør, ikke om bevegelsen skjer i et plan. En kloss på et skråplan beveger seg i et plan, men har bare én frihetsgrad.

Den dyreste enkeltfeilen i en 2D-kodedeloppgave: maks 2 av 5 poeng, uansett hvor korrekt resten av koden er.

Symptomet er lett å kjenne igjen: koden har bare y og vy, og beskriver dermed et loddrett kast i stedet for et skrått.

Fellen er nummer #6 i bokas feilregister. Motgiften er en mekanisk vane: tell tabellene før du skriver løkka. To romlige retninger krever fire tabeller, to akselerasjonslinjer og fire oppdateringslinjer.