raket.html

> restart;with(Units[Standard]):

Warning, the assigned name polar now has a global binding

Warning, these protected names have been redefined and unprotected: *, +, -, /, <, <=, <>, =, Im, Re, ^, abs, arccos, arccosh, arccot, arccoth, arccsc, arccsch, arcsec, arcsech, arcsin, arcsinh, arctan, arctanh, argument, ceil, collect, combine, conjugate, convert, cos, cosh, cot, coth, csc, csch, csgn, diff, eval, evalc, evalr, exp, expand, factor, floor, frac, int, ln, log, log10, max, min, normal, root, round, sec, sech, shake, signum, simplify, sin, sinh, sqrt, surd, tan, tanh, trunc, type, verify

hoogteRaket:=0; snelheidRaket:=0; inhoudRaket:=1500*Unit(centimeter^3); volumeWater:=600*Unit(centimeter^3); massaFles := 120*Unit(gram);

> dichtheidWater := 1000*Unit(kg/meter^3);

Oppervlakte van de opening waaruit het water uit de raket spuit

> oppervlakteStop := 4.9*Unit(centimeter^2);

> begindruk := 4*Unit(bar);

Een tijdstap van 1 milliseconde geeft bijvoorbeeld aan dat we elke duizendste seconde hoogte, snelheid, enz. zullen berekenen.

> tijdstap := 1*Unit(ms);

Dit zorgt ervoor dat we (te?) veel informatie zullen krijgen over wat er precies met onze raket gebeurt:

> printlevel := 2;

Laten we eens kijken wat er de eerste 10 milliseconden van de vlucht gebeurt.

>

beginVolumeWater := volumeWater; #onthouden voor later
beginVolumeLucht := inhoudRaket - beginVolumeWater;
tijd := 0;
while  tijd/Unit(ms) < 10  do
  if  volumeWater/Unit(meter^3) > 0  then
    volumeLucht := inhoudRaket - volumeWater;
    # de druk hangt af van het volume, en wel als volgt:
    druk := begindruk * (beginVolumeLucht/volumeLucht)^1.4;
    # kracht veroorzaakt door uitgestoten water
    krachtWater := 2 * druk * oppervlakteStop;
    # totale massa van fles en water dat er nog in zit
    totMassa := massaFles + volumeWater*dichtheidWater;
    volumeWaterVerloren :=
                  sqrt(2*druk/dichtheidWater)*oppervlakteStop*tijdstap;
  else
    # alle water is er al uit
    totMassa := massaFles;
    volumeWaterVerloren := 0;
  end if;
  zwaartekracht := totMassa * Unit(acceleration_of_free_fall);
  snelheidsverandering := (krachtWater-zwaartekracht)/totMassa * tijdstap;
  #### nieuwe gegevens een tijdstap later
  hoogteRaket := hoogteRaket + snelheidRaket*tijdstap;
  snelheidRaket := snelheidRaket + snelheidsverandering;
  volumeWater := volumeWater - volumeWaterVerloren;
  if  volumeWater/Unit(meter^3) < 0  then
    volumeWater := 0;   # anders komen we verderop in de problemen
  end if;
  tijd := tijd + tijdstap;
od;

Na 10 ms is de snelheid:

> snelheidRaket;

Of omgerekend

> convert(snelheidRaket,units,km/h);

Na iets meer dan 50 milliseconden (dit is 1/20e van een seconde!) is de raket al haar water kwijt.

> printlevel := 0; # een beetje minder uitvoer graag

hoogteRaket:=0; snelheidRaket:=0; inhoudRaket:=1500*Unit(centimeter^3); volumeWater:=600*Unit(centimeter^3); massaFles := 120*Unit(gram);

>

beginVolumeWater := volumeWater; #onthouden voor later
beginVolumeLucht := inhoudRaket - beginVolumeWater;
tijd := 0;
while  tijd/Unit(ms) < 52  do
  if  volumeWater/Unit(meter^3) > 0  then
    volumeLucht := inhoudRaket - volumeWater;
    # de druk hangt af van het volume, en wel als volgt:
    druk := begindruk * (beginVolumeLucht/volumeLucht)^1.4;
    # kracht veroorzaakt door uitgestoten water
    krachtWater := 2 * druk * oppervlakteStop;
    # totale massa van fles en water dat er nog in zit
    totMassa := massaFles + volumeWater*dichtheidWater;
    volumeWaterVerloren :=
                  sqrt(2*druk/dichtheidWater)*oppervlakteStop*tijdstap;
  else
    # alle water is er al uit
    totMassa := massaFles;
    volumeWaterVerloren := 0;
  end if;
  zwaartekracht := totMassa * Unit(acceleration_of_free_fall);
  snelheidsverandering := (krachtWater-zwaartekracht)/totMassa * tijdstap;
  #### nieuwe gegevens een tijdstap later
  hoogteRaket := hoogteRaket + snelheidRaket*tijdstap;
  snelheidRaket := snelheidRaket + snelheidsverandering;
  volumeWater := volumeWater - volumeWaterVerloren;
  if  volumeWater/Unit(meter^3) < 0  then
    volumeWater := 0;   # anders komen we verderop in de problemen
  end if;
  tijd := tijd + tijdstap;
od;

> hoogteRaket; convert(snelheidRaket,units,km/h); convert(volumeWater,units,cm^3);

Met een beetje prutsen kan je zien dat na 52 milliseconden er nog 6.36 cm^3 water in de raket zit en na 53 milliseconden is al het water al op.

Laten we nog een simulatie doen:

> hoogteRaket:=0; snelheidRaket:=0; inhoudRaket:=1500*Unit(centimeter^3); volumeWater:=600*Unit(centimeter^3); massaFles := 120*Unit(gram);

Maar we voegen een stukje code toe om grafieken te krijgen ... Ik heb ook een regeltje toegevoegd dat ervoor zorgt dat we in stappen van 10ms rekenen van zodra alle water uit de raket weg is, anders zou de berekening te lang duren (het duurt ook na die aanpassing nog lang genoeg).

>

tijdstap := 1*Unit(ms):
beginVolumeWater := volumeWater: #onthouden voor later
beginVolumeLucht := inhoudRaket - beginVolumeWater:
snelheidLijst := NULL:
hoogteLijst := NULL:
tijd := 0:
while  hoogteRaket/Unit(m) >= 0  do
  if  volumeWater/Unit(meter^3) > 0  then
    volumeLucht := inhoudRaket - volumeWater;
    # de druk hangt af van het volume, en wel als volgt:
    druk := begindruk * (beginVolumeLucht/volumeLucht)^1.4;
    # kracht veroorzaakt door uitgestoten water
    krachtWater := 2 * druk * oppervlakteStop;
    # totale massa van fles en water dat er nog in zit
    totMassa := massaFles + volumeWater*dichtheidWater;
      volumeWaterVerloren :=
                  sqrt(2*druk/dichtheidWater)*oppervlakteStop*tijdstap;
  else
    # alle water is er al uit
    totMassa := massaFles;
    krachtWater := 0;
    volumeWaterVerloren := 0;
  end if;
  zwaartekracht := totMassa * Unit(acceleration_of_free_fall);
  snelheidsverandering := (krachtWater-zwaartekracht)/totMassa * tijdstap;
  
  #### nieuwe gegevens een tijdstap later
  hoogteRaket := hoogteRaket + snelheidRaket*tijdstap;
  snelheidRaket := snelheidRaket + snelheidsverandering;
  volumeWater := volumeWater - volumeWaterVerloren;
  if  volumeWater/Unit(meter^3) < 0  then
    volumeWater := 0;   # anders komen we verderop in de problemen
    tijdstap := 10*Unit(ms);
  end if;

  #### nieuwe gegevens in de lijsten bijschrijven
  snelheidLijst := snelheidLijst , [tijd/Unit(ms),snelheidRaket/Unit(km/h)];
  hoogteLijst := hoogteLijst , [tijd/Unit(ms),hoogteRaket/Unit(m)];

  tijd := tijd + tijdstap;
od:

We kijken even naar de verzamelde gegevens.

> snelheidLijst;

Een grafiekje zegt veel meer dan tien bladzijden cijfertjes:

> plot([snelheidLijst],labels=["tijd (ms)","snelheid (km/u)"]);

> plot([hoogteLijst],labels=["tijd (ms)","hoogte (m)"]);

Hoe lang duurde de vlucht precies?

> tijd;

80 meter is wel nogal hoog; in de praktijk gaat de raket lang niet zo hoog. Wat zou er mis kunnen zijn aan onze simulatie?

> hoogteRaket:=0; snelheidRaket:=0; inhoudRaket:=1500*Unit(centimeter^3); volumeWater:=600*Unit(centimeter^3); massaFles := 120*Unit(gram); begindruk:=4*Unit(bar);

We houden nu ook rekening met de luchtweerstand. Die hangt af van de oppervlaktedoorsnede van de raket (dit is de doorsnede van de raket zoals gezien van bovenaf) en ook van de vorm van de raket. Die vorm wordt verwerkt door de weerstandcoefficient die 1 is voor slecht gestroomlijnde raketten maar 0,1 kan zijn voor zeer gestroomlijnde raketten. Het bepalen van die coefficient is niet bepaald eenvoudig dus we slaan er maar even naar en gokken op 0,5.

>

tijdstap := 1*Unit(ms):
weerstandCoeff := 0.5:
oppervlakteDoorsnede := 0.01*Unit(m^2):
beginVolumeWater := volumeWater: #onthouden voor later
beginVolumeLucht := inhoudRaket - beginVolumeWater:
snelheidLijst := NULL:
hoogteLijst := NULL:
tijd := 0:
while  hoogteRaket/Unit(m) >= 0  do
  if  volumeWater/Unit(meter^3) > 0  then
    volumeLucht := inhoudRaket - volumeWater;
    # de druk hangt af van het volume, en wel als volgt:
    druk := begindruk * (beginVolumeLucht/volumeLucht)^1.4;
    # kracht veroorzaakt door uitgestoten water
    krachtWater := 2 * druk * oppervlakteStop;
    # totale massa van fles en water dat er nog in zit
    totMassa := massaFles + volumeWater*dichtheidWater;
    volumeWaterVerloren :=
                  sqrt(2*druk/dichtheidWater)*oppervlakteStop*tijdstap;
  else
    # alle water is er al uit
    totMassa := massaFles;
    krachtWater := 0;
    volumeWaterVerloren := 0;
  end if;
  zwaartekracht := totMassa * Unit(acceleration_of_free_fall);
  # de 0.6 kg/m^3 op de volgende regel is de helft van de dichtheid van lucht
  luchtweerstand := -0.6*Unit(kg/m^3)*weerstandCoeff*oppervlakteDoorsnede
       * abs(snelheidRaket) * snelheidRaket;
  snelheidsverandering := (krachtWater-zwaartekracht+luchtweerstand)/totMassa * tijdstap;
  
  #### nieuwe gegevens een tijdstap later
  hoogteRaket := hoogteRaket + snelheidRaket*tijdstap;
  snelheidRaket := snelheidRaket + snelheidsverandering;
  volumeWater := volumeWater - volumeWaterVerloren;
  if  volumeWater/Unit(meter^3) < 0  then
    volumeWater := 0;   # anders komen we verderop in de problemen
    tijdstap := 10*Unit(ms);
    ## doe het hekje op de volgende regel weg om enkel de uitstootfase te volgen
    # break;
  end if;

  ## doe het hekje op de volgende regel weg om enkel het stijgen te volgen
  # if  snelheidRaket/Unit(meter/s) < 0  then break; end if;

  #### nieuwe gegevens in de lijsten bijschrijven
  snelheidLijst := snelheidLijst , [tijd/Unit(ms),snelheidRaket/Unit(km/h)];
  hoogteLijst := hoogteLijst , [tijd/Unit(ms),hoogteRaket/Unit(m)];

  tijd := tijd + tijdstap;
od;

> plot([snelheidLijst],labels=["tijd (ms)","snelheid (km/u)"]);

> plot([hoogteLijst],labels=["tijd (ms)","hoogte (m)"]);

Door de luchtweerstand is de hoogte een stuk lager en duurt de vlucht ook korter:

> tijd;

>