Tykillä ampuminen 2

Mallinnettaessa heittoliikettä, kuten esimerkiksi tykillä ampumista, keskeisinä vaikuttavina tekijöinä ovat painovoima sekä ilmanvastuksen aiheuttama nopeudelle vastakkaissuuntainen voima.

Painovoima on , missä on maan vetovoiman aiheuttama kiihtyvyys ja m heitettävän kappaleen massa. Ilmanvastusvoima on verrannollinen nopeuden neliöön, ja se voidaan mallintaa lausekkeella

missä b on ilmanvastuskerroin, v skalaarinen nopeus, nopeus vektorina ja nopeuden suuntainen yksikkövektori.

Yleinen liikeyhtälö saa tällöin muodon

missä on merkitty sekä .

Tarkastellaan tilannetta, jossa tykillä ammutaan poispäin liikkuvaan maaliin. Mikä on oikea ampumiskulma, kun 25 metrin korkeudella sijaitsevalta tykkilavetilta ammutaan 10 kg kranaatti nopeudella 1000 m/s kohti maalia, jonka etäisyys ampumishetkellä on tasan 5 km? Maali liikkuu 15 m/s tykistä poispäin veden pintaa pitkin, ja ilmanvastuskerroin kranaatille on .

Jaetaan yllä esitetty toisen kertaluvun vektorimuotoinen differentiaaliyhtälö kahteen eri osayhtälöön, z-suuntaan ja x-suuntaan, ja ratkaistaan yhtälöt. Laskujen aluksi on syytä hävittää mahdollisista aiemmista laskuista jääneet muuttujat.

> restart;

Määritellään differentiaaliyhtälöryhmä ja sen tuntemattomat muuttujat.

> ryhma:= m*diff(x(t), t$2)=-b*diff(x(t), t)*sqrt(diff(x(t), t)^2+diff(z(t), t)^2), m*diff(z(t), t$2)=-m*g-b*diff(z(t), t)*sqrt(diff(x(t), t)^2+diff(z(t), t)^2);

> muuttujat:= x(t), z(t);

Määritellään vakiot ja alkuehto. Ammuksen lähtökulma olkoon . Sisällytetään ehtoon sekä ammuksen alkunopeus jaettuna nopeuskomponentteihin että lähtökorkeus.

> g:= 9.81: b:= 10^(-4): m:= 10: v0:= 1000:

> alkuehto:= x(0)=0, z(0)=25, D(x)(0)=v0*cos(theta), D(z)(0)=v0*sin(theta);

Differentiaaliyhtälöryhmä ei ole ratkaistavissa alkeisfunktioiden avulla, joten ryhmää ei voi ratkaista Maple n avulla :n funktiona. Suljetussa muodossa olevan ratkaisun sijaan joudutaan turvautumaan arvaukseen. Arvataan jokin kulma (radiaaneissa), ja katsotaan, miten hyvin laukaus osui.

> theta:= 0.03;

> rtk1:= dsolve({ryhma, alkuehto}, {muuttujat}, numeric, startinit=true, output=listprocedure);

Poimitaan ratkaisusta x- ja z-suuntia kuvaavat komponentit listaksi ja piirretään kuvio.

> lentorata:= subs(rtk1, [x(t), z(t)]):

> plot([lentorata[], 0..10]);

Katsotaan miten lähelle laukaus osui. Huomaa, että ratkaisussa on huomioitu myös kohteen liike.

> lentoaika:= fsolve('lentorata[2](t)=0', t=0..10);

> osmaet:= lentorata[1](lentoaika)-(5000+15*lentoaika);

Ammus osuu veteen noin puolitoista kilometriä liian kauaksi. Lähtökulma oli siis liian suuri. Arvataan uudestaan tarkoituksena haarukoida oikea tulos arvausten väliin. Tämän jälkeen voidaan systemaattisesti hakea oikea kulma.

> theta:= 0.01;

> rtk2:= dsolve({ryhma,alkuehto}, {muuttujat}, numeric, startinit=true, output=listprocedure);

> lentorata:= subs(rtk2, [x(t), z(t)]):

> lentoaika:= fsolve('lentorata[2](t)=0', t=0..10);

> osumaet:= lentorata[1](lentoaika)-(5000+15*lentoaika);

Nyt laukaus jäi liian lyhyeksi. Saimme kuitenkin haarukoitua osumisvälin. Löytääksemme tarkan ampumiskulman, teemme organisoidun kokeilun.

Tehdään pieni for -silmukka, joka tulostaa osumapaikan 21 laskupisteessä 0.001 radiaanin välein.

> indeksi:= 0:
for theta from 0.01 by 0.001 to 0.03 do
rtk:= dsolve({ryhma, alkuehto}, {muuttujat}, numeric, startinit=true, output=listprocedure);
lentorata:= subs(rtk, [x(t), z(t)]);
lentoaika:= fsolve('lentorata[2](t)=0', t=0..10);
osumataul[indeksi]:= lentorata[1](lentoaika)-(5000+15*lentoaika);
kulmataul[indeksi]:= theta;
indeksi:= indeksi+1;
od:

kulmat:=[seq(kulmataul[k], k=0..20)]:
etaisyydet:= [seq(osumataul[k],k=0..20)];

Etsitään lähintä osumaa vastaava indeksi ja tämän avulla ampumiskulma ja osumaetäisyys.

> member(min(map(x->abs(x),etaisyydet)[]),etaisyydet,'ind'):
lahin:=[kulmat[ind],etaisyydet[ind]];

Sama uudestaan uudella kapeammalla haarukalla ja 0.0001 radiaanin askelvälillä.

> indeksi:= 0:
for theta from 0.020 by 0.0001 to 0.022 do
rtk:= dsolve({ryhma, alkuehto}, {muuttujat}, numeric, startinit=true, output=listprocedure);
lentorata:= subs(rtk, [x(t), z(t)]);
lentoaika:= fsolve('lentorata[2](t)=0', t=0..10);
osumataul[indeksi]:= lentorata[1](lentoaika)-(5000+15*lentoaika);
kulmataul[indeksi]:= theta;
indeksi:= indeksi+1;
od:

kulmat:=[seq(kulmataul[k], k=0..20)]:
etaisyydet:= [seq(osumataul[k],k=0..20)];

> member(min(map(x->abs(x),etaisyydet)[]),etaisyydet,'ind'):
lahin:=[kulmat[ind],etaisyydet[ind]];

Osuma on siis runsaan kuuden metrin päässä ja vastaava ampumiskulma on asteissa

> convert(lahin[1], degrees):
evalf(%);

Tehtäviä

Piirrä lähimmäksi osuvan ammuksen lentorata. Määritä radan lakikorkeus ja ammuksen lentoaika.

Suurenna ilmanvastuskerrointa, ja tutki, miten se vaikuttaa radan muotoon.