Tykillä ampuminen 2

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

Painovoima on \(G\&_ = m\ g\&_\), missä \(g\&_\) on maan vetovoiman aiheuttama kiihtyvyys ja m heitettävän kappaleen massa. Ilmanvastusvoima on verrannollinen nopeuden neliöön, ja se voidaan mallintaa lausekkeella

\(\(\(F\&_ = \(\(-b\)\ \(v\^2\) 
          v\&^\  = \ \(-b\)\ v\ v\&_\)\)\(,\)\)\)

missä b on ilmanvastuskerroin, v skalaarinen nopeus, \(v\&_\) nopeus vektorina ja \(v\&^\) nopeuden suuntainen yksikkövektori.

Yleinen liikeyhtälö \(F\&_ = m\ a\&_\) saa tällöin muodon

\(\(\(m \(\( d\^2\) r\&_\)\/\(d\ t\^2\) = 
      m  g\&_\  - 
          b\  | \(d\ r\&_\)\/\(d\ t\) | \ \(d\ r\&_\)\/\(d\ t\)\)\(\ \)\(,\)\)\)

missä on merkitty \(a\&_\)=\(\(\(d\^2\) r\&_\)\/\(d\ t\^2\)\) sekä \(v\&_\)=\(\(d  r\&_\)\/\(d\ t\)\).

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  b=(10^(-4)).

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.

Remove["Global`*"]

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

ryhma =  {m x''[t] == -b x'[
              t] Sqrt[x'[t]^2 + z'[t]^2],  
        m z''[t]  == -m g - 
            b z'[
                t] Sqrt[x'[t]^2 + z'[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,  x'[0] == v0 Cos[θ],  z'[0] == 
          v0 Sin[θ]}

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

rtk1 =  
      NDSolve[Flatten[{ryhma, (alkuehto /. θ -> 0.03)}], 
        muuttujat, {t, 0, 10}]

Muokataan tulos käyttökelpoiseen muotoon ja piirretään kuvio.

lentorata = muuttujat/.First[
          rtk1]

ParametricPlot[lentorata, {t, 0, 10}]

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

osumaet = 
      lentorata[[1]] - (5000 + 15  t) /. 
        FindRoot[lentorata[[2]] == 0, {t, 0, 10}]

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.

rtk2 =  
        NDSolve[Flatten[{ryhma, (alkuehto /. θ -> 0.01)}], 
          muuttujat, {t, 0, 10}];
lentorata = 
        muuttujat /. First[rtk2] ;
osumaet = 
      lentorata[[1]] - (5000 + 15  t) /. 
        FindRoot[lentorata[[2]] == 0, {t, 0, 10}]

Nyt laukaus jäi liian lyhyeksi. Oikea lähtökulma on kuitenkin saatu haarukoiduksi. Tarkan ampumiskulman löytämiseksi tehdään organisoitu kokeilu, johon otetaan aluksi 21 laskupistettä 0.001 radiaanin välein.

Lasketaan kulma ja sen avulla lasketut yhtälöt omaan taulukkoonsa ja osumaetäisyydet omaansa.

lentoradat = 
        Table[Flatten[{θ, 
              muuttujat /. 
                First[NDSolve[Flatten[{ryhma, alkuehto}], 
                    muuttujat, {t, 0, 10}]]}], {θ, 0.01, 0.03, 
            0.001}];

etaisyydet = 
      Table[lentoradat[[i, 2]] - (5000 + 15  t) /. 
          FindRoot[lentoradat[[i, 3]] == 0, {t, 0, 10}], {i, 1, 
          21}]

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

indeksi =  
        Select[Range[21], 
          Abs[etaisyydet[[#]]]  ==  
              Min[Abs[etaisyydet]]  &];
lahin =  {lentoradat[[indeksi, 1]], 
        etaisyydet[[indeksi]]}

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

lentoradat = 
        Table[Flatten[{θ, 
              muuttujat /. 
                First[NDSolve[Flatten[{ryhma, alkuehto}], 
                    muuttujat, {t, 0, 10}]]}], {θ, 0.020, 
            0.022, 0.0001}];

etaisyydet = 
      Table[lentoradat[[i, 2]] - (5000 + 15  t) /. 
          FindRoot[lentoradat[[i, 3]] == 0, {t, 0, 10}], {i, 1, 
          21}]

indeksi =  
        Select[Range[20], 
          Abs[etaisyydet[[#]]]  ==  
              Min[Abs[etaisyydet]]  &];
lahin =  {lentoradat[[indeksi, 1]], 
        etaisyydet[[indeksi]]}

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

lahin[[1, 1]]/Degree

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.

JP & SKK 03.05.2001