laskuv.mws |
Heittoliikkeeseen kuuluu keskeisinä tekijöinä maan vetovoima sekä ilmanvastuksen aiheuttama nopeudelle vastakkaissuuntainen voima. Näitä voimia sovelletaan myös mallinnettaessa laskuvarjolla hyppäämistä.
Tarkastellaan ensin vapaata putoamista. Putoajaan vaikuttava voima aiheuttaa kiihtyvyyden, joka saa aikaan nopeuden. Yleisessä muodossa voimalle pätee
,
missä m on tarkasteltavan kohteen massa ja kiihtyvyys. Hyppääjään vaikuttavat voimat ovat siis painovoima , missä on maan vetovoiman aiheuttama kiihtyvyys, sekä ilmanvastus, jota arvioidaan yhtälöllä
Tässä b on ilmanvastuskerroin, v skalaarinen nopeus, nopeus vektorina ja nopeuden suuntainen yksikkövektori. Liikeyhtälö saa tällöin muodon
missä on merkitty ja .
Tämä on heittoliikkeen differentiaaliyhtälö yleisessä muodossa.
Tehdään tarvittavat laskut laskuvarjohyppääjän lentoradan ja -nopeuden selvittämiseksi. Edellä esitetty ilmanvastusmalli on vain approksimaatio todelliselle ilmanvastukselle. On myös huomattava, että mikäli hyppääjä putoaa ilman varjoa, ilmanvastustekijä on aivan eri suuruinen kuin varjon kanssa.
Laskujen aluksi on syytä hävittää mahdollisista aiemmista laskuista jääneet muuttujat.
> restart;
Määritellään differentiaaliyhtälö normaaliryhmänä, jossa on erikseen x ja y-komponentit.
> ryhma:= diff(x(t), t)=u(t), m*diff(u(t), t)=-b*sqrt(u(t)^2+v(t)^2)*u(t), diff(y(t), t)=v(t), m*diff(v(t), t)=-m*g-b*sqrt(u(t)^2+v(t)^2)*v(t);
Probleeman tuntemattomat funktiot ovat vaakasijaintia kuvaava x( t ) sekä korkeus y( t ) ja vastaavat nopeudet u( t ) ja v( t ).
> tuntemattomat:= x(t), y(t), u(t), v(t);
Määritetään parametrit ja valitaan alkuehdot lentokoneesta hypättäessä. Oletetaan hitaasti lentävä koulutuslentokone.
> m:= 70: b:= 30: g:= 9.81:
> alkuehto:= x(0)=0, y(0)=1000, u(0)=20, v(0)=0;
Ratkaistaan differentiaaliyhtälöryhmä. Käytetään startinit = true -parametria, koska myöhemmin käytettävä fsolve -komento vaatii sitä.
> rtk:= dsolve({ryhma,alkuehto}, {tuntemattomat}, type=numeric, output=listprocedure, startinit=true);
Poimitaan ratkaisusta hyppääjän lentorataa ja nopeutta kuvaavat komponentit listaksi.
>
rata:= subs(rtk, [x(t), y(t)]):
nopeus:= subs(rtk, [u(t), v(t)]):
Piirretään kuva laskuvarjohyppääjän lentoradasta. Huomaa koordinaattiakselien eri skaalat!
> plot([rata[], 0..250], view=[0..10, 950..1000]);
> fsolve('rata[2](t)=0');
Lentorata ei muodostunut kovin pitkäksi ja laskeutuminen kilometrin korkeudelta tapahtui n. 3.5 minuutissa. Tarkastellaan vielä hyppääjän nopeutta.
> plot(nopeus, 0..250, view=[0..5, -6..20]);
> array([seq(nopeus(k), k=1..10)]);
Vaakasuuntainen nopeus lähes häviää ja putoamisnopeus asettuu vajaaseen 4.8 m/s nopeuteen jo kahden ensimmäisen sekunnin aikana. Hyppääjän putoamisvauhti maahan on siis noin 17 km/h, joka on suurenpuoleinen jaloilla vastaanotettavaksi.
Tarkastellaan seuraavaksi tilannetta, jossa hyppääjä ei avaa varjoaan välittömästi. Tehdään vastaava tarkastelu kuin edellä, paitsi että ilmanvastuskerroin
b
muuttuu hyppääjän avatessa varjon. Oletetaan, että ilman varjoa hyppäävän ihmisen maksiminopeus (vaikuttavien voimien summa on 0) eli ns. terminaalinopeus on runsaat 200 km/h eli n. 60 m/s. Tällöin ilmanvastuskerroin voidaan ratkaista yhtälöstä
= 0, jolloin
.
> b0:= 70*9.81/60^2;
Jotta voimme muuttaa ilmanvastuskertoimen arvon = 0.19075 varjon auetessa arvoon b = 30, tarvitaan funktio, joka mallintaa muutoksen. Tällaiseksi sopii
> f:= x->piecewise(x>-infinity, 1-exp(-exp(2*x-0.5)));
Funktio on määritelty paloittain, koska Maple versio 6 ei suoriutunut myöhemmin esiintyvän differentiaaliyhtälöryhmän ratkaisemisesta käytettäessä perinteisesti määriteltyä funktiota.
Funktion kuvaaja näyttää seuraavalta:
> plot(f(x), x=-10..10, view=[-10..10, -0.5..1.5]);
Määritetään putoamisyhtälöt uudestaan. Ilmanvastuskerroin määräytyy nyt funktiosta
> bb:= t->b0+(b-b0)*f(t);
ja uusi yhtälöryhmä on
> ryhma2:= diff(x(t), t)=u(t), m*diff(u(t), t)=-bb(t-20)*sqrt(u(t)^2+v(t)^2)*u(t), diff(y(t), t)=v(t), m*diff(v(t), t)=-m*g-bb(t-20)*sqrt(u(t)^2+v(t)^2)*v(t);
Varjon on oletettu aukeavan 20 sekuntia hyppäämishetken jälkeen.
Ratkaistaan yhtälöryhmä numeerisesti. Kyseessä oleva yhtälöryhmä on hankala ratkaistava numeerisesti ja oletusarvoilla dsolve -komento suoranaisesti laskee väärin. Startinit=true -määreellä laskenta perustuu aina annettuihin alkuarvoihin. Normaalisti laskenta perustuu jo laskettuihin välituloksiin. Startinit=true siis hidastaa laskentaa melkoisesti, koska laskenta aloitetaan aina alkuarvoista, mutta tulos on kuitenkin tarkempi.
> rtk2:= dsolve({ryhma2,alkuehto}, {tuntemattomat}, output=listprocedure, type=numeric, startinit=true);
Poimitaan ratkaisusta hyppääjän lentorataa ja nopeutta kuvaavat komponentit listaksi.
>
rata2:= subs(rtk2, [x(t), y(t)]):
nopeus2:= subs(rtk2, [u(t), v(t)]):
Piirretään kuvaaja hyppääjän lentoradasta, määritetään alastulohetki ja alastulonopeus:
> plot([rata2[], 0..50], view=[0..200, 0..1000]);
> alas:=fsolve('rata2[2](t)=0',t,0..50);
> nopeus2(alas);
Tarkastellaan vielä putoamisnopeuksia kokonaistilanteen hahmottamiseksi. Putoamisnopeudet ajan funktiona ovat
> with(plots):
Warning, the name changecoords has been redefined
> display(plot(nopeus2, 0..50, view=[0..50, -60..20]), plot([[17, y, y=-60..20], [20, y, y=-60..20]], color=black));
> rata2(17);
Putoaminen tapahtuu varsin eri tavalla kuin avoimella varjolla. Vaikka varjon avautumishetkeksi asetettiin 20 s, avautuminen on alkanut askelfunktion mukaisesti muuttamaan ilmanvastuskerrointa radikaalisti ja siten jarruttamaan putoamista jo aikaisemmin. Nopeus alkaa vähetä noin 250 metrin korkeudessa, kun t = 17 s.
Piirretään kuvaaja kiihtyvyydelle. Koska ratkaisuproseduuria ei voi derivoida, lasketaan derivaatalle arvioita käyttämällä erotusosamäärää.
Poimitaan aluksi ratkaisuproseduurista nopeuden pystysuuntainen komponentti:
> numv:= subs(rtk2, v(t));
Muodostetaan aikavälistä 0 - 50 lista 0.5 sekuntin välein.
> aika:= seq(t/2, t=0..100):
Lasketaan nopeuden erotusosamäärä aika-listan pisteissä.
> derivpisteet:= seq((numv(t/2+0.00001)-numv(t/2))/0.00001, t=0..100):
Piirretään kuvaaja derivaatasta ajan funktiona käyttäen edellä laskettuja erotusosamääriä.
> display(plot(zip((x, y)->[x, y], [aika], [derivpisteet])), plot([19, t, t=-10..30], color=black));
> (numv(19+0.00001)-numv(19))/0.00001;
Voimakkain jarrutuskiihtyvyys (n. 28 ) tapahtui kohdassa t = 19 s. Putoamisen kokonaisajasta avatun varjon kanssa kuluu yhä yli puolet, vaikka matka on noin viidennes.
Milloin hyppääjän on viimeistään vedettävä varjon avausnarusta, kun oletetaan, että tämän on tapahduttava 3 sekuntia ennen askelfunktiossa käytettyä aikaa? Suurin mahdollinen nopeus maahan pudotessa ilman loukkaantumista on 10 m/s. Mikä on tällöin putoamisen kokonaisaika?