De vlucht van een pijl zonder en met luchtweerstand

From Pijl en Boog
Jump to navigationJump to search

De vlucht van een pijl zonder luchtweerstand

  • Het traject van een pijl is, bij afwezigheid van luchtweerstand, enkel afhankelijk van de snelheid bij vertrek, de hoek waarin de pijl afgeschoten wordt, de hoogte van waarop de pijl vertrekt, de valversnelling en de hoogte van het eindpunt.
  • Om de formule te bepalen hebben we dus nodig:
    • snelheid
    • valversnelling
    • de hoek waaronder de pijl afgeschoten wordt (in graden)
    • begin hoogte (in m)
    • eind hoogte (in m)
    • : horizontale beginpositie
    • = de afgelegde horizontale afstand
    • : verticale beginpositie
    • = de afgelegde verticale afstand (=hoogte)
  • De positie van de pijl wordt bepaald door de formule (uit "Doing Physics with Scientific Notebook: A Problem Solving Approach, Joseph Gallant, 2012"):


  • de maximale afstand, R, die een pijl kan afleggen (zonder luchtweerstand) wordt gegeven door de formule

  • Wanneer de begin en eind hoogte dezelfde is, kan formula vereenvoud worden tot:

    • Deze afstand is maximaal voor een hoek van 45°. Dan wordt gelijk aan 1 en is de formule

  • De tijd () dat de pijl onderweg is wordt gegeven door de formule

  • De maximale hoogte () dat een pijl kan bereiken wordt gegeven door de formule:

  • Hieronder staan de trajecten voor een pijl afgeschoten aan 200 fps (links) en 270 fps (rechts).
    • De vorm van de trajecten is dezelfde voor de 200 fps en de 270 fps pijl, enkel de hoogte en afstand is verschillend.


  • Indien er geen luchtweerstand is, zou je zelfs met een moderne compound boog van 70# gemakkelijk over de Eifeltoren heen kunnen schieten!
    • Vandaar altijd zorgen voor een goede pijlenopvang achter het doel, nooit de boog naar boven aanspannen, enz.
  • Vermits een model zonder luchtweerstand niet realistisch is, dienen we de berekeningen te herhalen met luchtweerstand. (zie verder)

Vlucht van de pijl rekening houdend met de luchtweerstand.

  • De luchtweerstand voor een object in beweging wordt weergegeven door de volgende formule
  • met:
    • de dichtheid van het medium waarin het voorwerp zich voortbeweegt, hier de lucht.
    • de snelheid van het voorwerp in m/s
    • de drag coefficiënt (weerstandscefficiënt)
    • de oppervlakte van de doorsnede loodrecht op de bewegingsrichting
    • de kracht die uitgeoefend wordt op het voorwerp tijdens de beweging in het medium (lucht)
  • Hierboven hielden we enkel rekening met de valversnelling in de y-component, moeten we hier rekening houden met de luchtweerstand voor zowel de x- als de y-component van de snelheid.
    • beginsnelheid
    • de hoek waaronder de pijl wordt afgeschoten.
    • de x-component van de beginsnelheid
    • de x-component van de beginsnelheid
  • We kunnen dus op zowel de x- als de y-component van de snelheid, per klein tijdseenheid dat het voorwerp bewogen heeft, de nodige aanpassingen doen. Dit is de methode van Euler en varianten zoals Runge-Kutta en dit zijn numerieke methoden om zulk probleem op te lossen.
    • Per klein tijdsinterval kunnen we de aangepaste snelheid berekenen.
  • Vermits wil dit zeggen dat
    • Dus voor kunnen we dus stellen dat de (negatieve) versnelling door de luchtweerstand gelijk is aan:
    • Deze versnelling (=vertraging want werkt in tegengestelde richting aan de snelheid van het voorwerp) kunnen we dan gebruiken om de snelheid
      • voor de x-component van de snelheid
      • voor de y-component van de snelheid
        • hier werkt ook de valversnelling mee
      • Dus de nieuwe snelheid op verder in de tijd wordt dus berekend op basis van de voorgaande snelheid + de kleine aanpassing. Door dit te doen vanaf de beginsnelheid en dit te herhalen over het gehele traject bekomen we dus de baan van de pijl met luchtweerstand.
    • Een voorwerp dat naar beneden valt zal versnellen totdat het aan zijn eindsnelheid (terminal velocity) komt.. vanaf dat moment is de snelheid van het vallend voorwerp constant en versnelt het niet meer.
    • De Eindsnelheid/Terminal velocity wordt weergegeven door de formule:
    • In het computer programma dat we gebruiken om de baan van de pijl met luchtweerstand te berekenen, moeten we ook rekening houden dat indien de snelheid bij het terug naar beneden vallen groter zou worden dat de Eindsnelheid, om deze dus te beperken tot maximaal de Eindsnelheid/terminal velocity.
  • Wanneer we dit toepassen voor hoeken van 89, 75, 45 en 30 graden resulteert dit in volgende grafiek:
  • Vlucht van de pijl met en zonder luchtweerstand, uitgezet in de tijd.
  • Hoe zit het nu met hoger schieten dan de Eifeltoren... lukt dat dan nog met een boog en hoe snel moet je pijl dan wel niet zijn?
    • Spoiler: een 700 grain pijl aan 305 fps afgeschoten onder een hoek van 86° zou er nog net (332m) over vliegen en aan de andere kant van de voet van de Eifeltoren landen. Warning: do not try this at home.. or in Paris!

Nota over de drag coefficiënt

    • De -waarde van een pijl is sterk afhankelijk van hoe de pijl is opgebouwd. Een flu-flu pijl zal uiteraard veel sneller remmen en dus een veel grotere -waarde hebben dan een dan een gewone pijl en kleine veren die perfect parallel aan de schacht gekleefd zijn.
    • Een hoge spine-waarde, en dus een slappe spine, zal de pijl doen "paradoxen" over z'n traject dus ook de -waarde laten toenemen.
    • In een artikel werd de -waarde van de pijl berekend door de pijl magnetisch op te hangen in een windtunnel. Dit levert een zeer accurate meting op van de -waarde van de pijl. In het artikel vonden ze -waarden van 2.6 voor een pijl met spin-wings en een -waarde van 1.5 voor een pijl met GasPro vanes.
    • Een naakte schacht heeft een -waarde van zoals getoond in deze video dicht bij de waarde van 0.94 uit dit artikel.
    • De -waarde die we hier gebruiken van is afkomstig van een video van een persoon die een aantal verschillende pijlen onder een hoek van 45° afschiet.
    • Ter vergelijking, een vrachtauto heel een -waarde van 0.6, dus minder dan de helft van sommige pijlen.

Onder welke hoek kunnen we het verste schieten

Met een pijl

Waar voor het model zonder luchtweerstand je het verste kan schieten door onder een hoek van 45° te schieten, is dit voor het model met luchtweerstand niet zo.

Zo is het bv dat voor een pijl van 550 grain @ 270 fps je het verste kan schieten door onder een hoek van 42° te schieten.

Je kan dit testen door in de formule telkens een andere hoek in te geven. Er bestaan uiteraard wiskundige manieren voor om dit optimum te berekenen.

Met een zware M114 155 mm howitzer kogel

Voor een M114 155 mm howitzer kogel schiet je het verste onder een hoek van 30.5°. Dit in de veronderstelling dat de drag coefficiënt 0.295 is.

Python code

  • Syntaxhighlighting lijkt niet te werken, sorry.

<syntaxhighlight lang="python" line> import numpy as np print ("Hello archery world!") ‎</syntaxhighlight>

imports


import numpy as np
import matplotlib.pyplot as plt

helper functies

# Define helper functions
# KE
def kinetic_energy(m, v):
    return 0.5*(m*0.0647989/1000)*(v*0.3048)*(v*0.3048)
# momentum
def momentum(m, v):
    return (m*0.0647989/1000)*(v*0.3048)

def grain2gram(x):
    return x * 0.0647989

def grain2kg(x):
    return grain2gram(x)/1000

def gram2grain(x):
    return x / 0.0647989

def kg2grain(x):
    return x*15432.3584

def fps2ms(x):
    return x * 0.3048

def ms2fps(x):
    return x / 0.3048

def inch2meter(x):
    return(x*0.0254)
def meter2inch(x):
    return(x/0.0254)

def lbs2Newton(x):
    return(x*4.4482216)
def Newton2lbs(x):
    return(x/4.4482216)

Functies die het traject van de pijl met en zonder luchtweerstand berekenen

#####################
# PIJLENBOOG MASTER #
#####################

import numpy as np
import matplotlib.pyplot as plt

def simuleer_traject(x0, y0, v0, theta, dt, num_timesteps, diameter, C_factor, air_density,  arrow_mass, AR, stop_height):
    """
    Simuleert de baan van een pijl onder invloed van luchtweerstand.
    
    Parameters:
    x0 (float): begin x-coördinaat van de pijl
    y0 (float): begin y-coördinaat van de pijl: begin hoogte
    v0 (float): begin snelheid van de pijl
    theta (float): initiële lanceerhoek van de pijl (in radialen)
    dt (float): tijd stapgrootte
    num_timesteps (int): aantal tijdstappen om te simuleren
    air_density (float): luchtdichtheid
    AR: (boolean) simuleer luchtweerstand (True/False)
    stop_height (float): stop hoogte
     
    Returns:
    xs (list of floats): x-coordinates of the arrow at each timestep
    ys (list of floats): y-coordinates of the arrow at each timestep
    timevalues (list of floats): tijdwaarden in seconden verlopen na lancering
    """
    # valversnelling
    g=9.81
    
    # bereken de sectional diameter van de pijl
    A=np.pi*(diameter/2)**2
    
    # oppervlakte van 3 blazer vanes
    # Dikte = 0.4 mm, hoogte 13.5 mm
    # (0.4/1000)*(13.5/1000) m²
    #
    verenoppervlak=3*(0.4/1000)*(13.5/1000)
    # Totaal sectional oppervlak van de pijl.
    A=A+verenoppervlak
  
    # Instellen van de beginwaarden
    xs = [x0]
    ys = [y0]
    vs = [v0]
    # de x- en y-component van de snelheid
    vx = v0 * np.cos(theta)
    vy = v0 * np.sin(theta)
    v=v0
    vxs = [vx] # list met de x-component van de snelheid (horizontaal)
    vys = [vy] # List met de y-component van de snelheid (horizontaal)
    # de lijst met versnellingen in de x- en y-richting
    axs=[]
    ays=[]
    # Indien luchtweerstand in rekening moet gebracht worden
    if AR==True:
        # Formule van luchtweerstand
        FA= 0.5*air_density*A*C_factor*v**2
        # De x- en y-component van de luchtweerstand. Die heeft dezelfde verhouding als de x en y component van de snelheid. Negatief omdat die tegengesteld is aan de snelheid.
        FAx= -FA*vx/v
        FAy= -FA*vy/v
        # F=m.a. Dus de versnelling is de kracht gedeeld door de massa. 
        ax = FAx/arrow_mass
        ay = FAy/arrow_mass-g
        axs.append(ax)
        ays.append(ay)
        
    else:
        # hier is geen luchtweerstand dus telt enkel de valversnelling.
        ax=0
        ay=-g
        axs.append(ax)
        ays.append(ay)

    # De delta vx, delta vy, delta x en delta y
    dvxs=[]
    dvx=ax*dt
    dvxs.append(dvx)
    
    dvys=[]
    dvy=ay*dt
    dvys.append(dvy)
    
    dxs=[]
    dx=(vx + dvx/2)*dt # dvx/2 is de gemiddelde x-snelheid in het dt interval
    dxs.append( dx )
    
    dys=[]
    dy=(vy + dvy/2)*dt # dvy/2 is de gemiddelde y-snelheid in het dt interval
    dys.append( dy )
    
    v_term=0
    if AR==True:
        # Het is niet heel waarschijnlijk dat we voor een pijl de Eindsnelheid/terminal velocity zullen bereiken, maar op deze manier is de functie ook bruikbaar voor andere zaken.
        v_term=np.sqrt(2*arrow_mass*g/(air_density*A*C_factor))
    
    timevalues=[0]
    
    for i in range(num_timesteps):
        # Update van de nieuwe positie en snelheid
        xs.append(xs[-1] + dxs[-1])
        ys.append(ys[-1] + dys[-1])
        
        vx=vxs[-1]+dvxs[-1]
        vy=vys[-1]+dvys[-1]
        vxs.append(vx)
        vys.append(vy)

        v=np.sqrt(vx**2+vy**2)
        vs.append(v)
        
        # Berekening van de luchtweerstand en de versnelling in x en y richting.
        if AR==True:
            FA= 0.5*air_density*A*C_factor*v**2
            FAx= -FA*vx/v
            FAy= -FA*vy/v
            ax = FAx/arrow_mass
            ay = FAy/arrow_mass-g
            axs.append(ax)
            ays.append(ay)
        else:
            ax=0
            ay=-g
            axs.append(ax)
            ays.append(ay)

        dvx=ax*dt
        dvy=ay*dt
        dvxs.append(dvx)
        dvys.append(dvy)
                
        dx=(vx +  dvx/2)*dt # dvx/2 is de gemiddelde x-snelheid in het dt interval
        dy=(vy +  dvy/2)*dt # dvy/2 is de gemiddelde y-snelheid in het dt interval
        
        dxs.append(dx)
        dys.append(dy)
        
        
        # wanneer de neerwaardse snelheid groter wordt dan de terminal velocity, maak het gelijk aan de terminal velocity.
        if (vy < -1*v_term and AR==True):
            vy= -1*v_term

        timevalues.append(timevalues[-1]+dt)
# We negeren Stokes nog even.. enkel bij trage snelheden...misschien later toevoegen
#         if (np.sqrt(vx**2+vy**2)*diameter)<0.015:
#             print("Stokes drag should be used!")
        if ((ys[-1] ) <= stop_height and vy<0):  # stop simulation when vy becomes smaller than or equal to zero.. let arrow stop at ground
            break
    return xs, ys, timevalues




# De functies hieronder krijgen de hoek in radialen en houden geen rekening met luchtweerstand.
def y_positie(hoek,snelheid,start_hoogte, eind_hoogte, start_x, eind_x):
    g=9.81
    theta=hoek
    return ( start_hoogte + np.tan(theta)*(eind_x-start_x) - g*(eind_x-start_x)**2/(2*(snelheid*np.cos(theta))**2  )  )

def max_distance(hoek,snelheid,start_hoogte,eind_hoogte):
    g=9.81
    theta=hoek
    return ( snelheid**2*np.sin(2*theta)/(2*g) * ( 1 + np.sqrt(1 - 2*g*(eind_hoogte-start_hoogte)/( (snelheid*np.sin(theta))**2)  )  )  )

Eigenlijk programma om de plot aan te maken.

#######################################################

# Parameters over de lanceerhoek, gewicht van de pijl, drag coefficiënt, diameter, luchtweerstand en snelheid van de pijl.
# Door deze parameters in een lijst te zetten kunnen we gemakkelijk verschillende grafieken maken door deze parameters aan te passen.

hoeken=[89,75,45,30] # hoek in graden
hh=list(np.array(hoeken)*np.pi/180) # omzetten naar radialen. Kan in formule via truuk met numpy array ;-)
mm=[grain2kg(550)]*len(hoeken) # gewichten in grains, we herhalen de waarde evenveel keer als dat er hoeken zijn. Indien nodig kunnen we ook verschillende gewichten ingeven.
# CC=[1.42]*len(hoeken) # drag coefficiënt
CC=[0.94]*len(hoeken) # drag coefficiënt. We gebruiken 0.94. Andere waarden zijn mogelijk in functie van de gekozen pijl.
dd=[inch2meter(0.280)]*len(hoeken) # diameter
rr=[1.3]*len(hoeken) #kg/m^3 Air density
ss=[fps2ms(270)]*len(hoeken)

# Verschillende waarden voor Cd getest (drag coefficient)
# C=0.295 # Bullet https://en.wikipedia.org/wiki/Drag_coefficient
# C=1.65 # https://www.sciencedirect.com/science/article/pii/S1877705813010680#:~:text=The%20arrow%20launched%20by%20the,induces%20the%20boundary%20layer%20transition.
# C=0.9 # based on video, no fletching, https://www.youtube.com/watch?v=1L8r4hAgVsw
# C=1.42 # experimental by someone in a field
# The shape of the projectile influences its coefficient of drag, which determines how strongly air resistance affects an object's flight path. The higher the coefficient, the larger the effect. 
# For shapes, you can choose "none," which has zero drag; "wing shape," such as the wing of an airplane with a very low coefficient of 0.2; a "sphere," which has a coefficient of 0.5; 
# and a "flat plate," which is anything with a face of a large, flat surface area, with a coefficient of 1.28.
# Uiteindelijk een waarde van 1.42 genomen.. dit moet experimenteel bevestigd worden. De Cd uit het artikel is voor een heel slappe pijl die veel paradox heeft, een stijve pijl uit een compound boog zal een lagere Cd hebben. 
# Spinwing vs Gaspro veren gaf ook een heel groot verschil.
# Een andere persoon berekend Cd aan de hand van de ratio met G1... Dit zou voor zo'n pijl Cd=4.1 opleveren, wat me veel lijkt, maar niet onmogelijk.
# https://sites.google.com/site/technicalarchery/technical-discussions-1/drag-coefficients-of-bullets-arrows-and-spears

# Stel de initiele voorwaarden in
g=9.81 #N/kg
# Start positie en gewenste hoogte om te stoppen.
launch_x = 0
launch_y = 0
stop_y=0

# Hoe dient de grafiek gemaakt te worden. Taal en doelpubliek en dient er een logo gebruikt te worden.
LANG='NL'

dt = 0.015  # Tijdsinterval 
num_timesteps = 10000  # simulatie voor 1000 tijdsintervallen

# De lijst met gebruikte hoeken voor gebruik in de titel van de grafiek.
hoeken_str=''
for i in np.arange(0,len(hoeken)-2):
    hoeken_str=hoeken_str+ str(hoeken[i]) + '°, '
hoeken_str=hoeken_str+ str(hoeken[len(hoeken)-2])+'° '
if LANG=="EN":
    hoeken_str=hoeken_str+ 'and ' + str(hoeken[len(hoeken)-1])+'°'
else:
    hoeken_str=hoeken_str+ 'en ' + str(hoeken[len(hoeken)-1])+'°'

    
# Klaarmaken van de tekst voor titel en assen.
if LANG=='EN':
    my_x_axis_label='Distance (m)'
    my_y_axis_label='Height (m)'
    my_title='Trajectory of a {:.0f} grain arrow for '.format(kg2grain(mm[0]))+hoeken_str+' @ {:.0f} fps with and without air resistance ($C_d={:.2f}$)\nstarting at a height of {:.1f} m and arriving at a height of {:.1f} m with 3 blazer vanes'.format(ms2fps(ss[0]),CC[0],launch_y,stop_y)
elif LANG=='NL':
    my_x_axis_label='Afstand (m)'
    my_y_axis_label='Hoogte (m)'
    my_title='Vlucht van een {:.0f} grain pijl voor '.format(kg2grain(mm[0]))+hoeken_str+' @ {:.0f} fps met en zonder luchtweerstand ($C_d={:.2f}$)\nvertrekkend van {:.1f} m hoogte en eindigend op {:.1f} m hoogte met 3 blazer vanes'.format(ms2fps(ss[0]),CC[0],launch_y,stop_y)
else:
    my_x_axis_label='X'
    my_y_axis_label='Y'
    my_title='Title'
   
fig, ax = plt.subplots(figsize=(20, 10))
ax.set_title(my_title,fontsize=20)
ax.set_xlabel(my_x_axis_label,fontsize=16)
ax.set_ylabel(my_y_axis_label,fontsize=16)

# Stel de gewenste parameters van de grafiek in.
plt.xticks(np.arange(0, 860, 20) )
ymin=-5
ymax=375
ax.set_ylim(ymin, ymax)
ax.set_xlim(-5, 700)

# De volgorde van de gebruikte kleuren.
my_colors=[u'#1f77b4', u'#ff7f0e', u'#2ca02c', u'#d62728', u'#9467bd', u'#8c564b', u'#e377c2', u'#7f7f7f', u'#bcbd22', u'#17becf']
# Hoe dient de tekstbox er uit te zien.
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

# Simulatie traject voor elke opgegeven hoek en bijhorende parameters
for index,hoek in enumerate(hh):
    my_x, my_y, my_t= simuleer_traject(launch_x, launch_y, ss[index], hoek, dt, num_timesteps, dd[index] , CC[index], rr[index], mm[index],True,stop_y)

    plt.plot(my_x, my_y,color=my_colors[index])
    # Maak een test string met de hoek in graden en het punt (x,y) van de maximum.
    textstr = '{:.0f}°: max={:.0f}, {:.0f}'.format(hoeken[index], my_x[my_y.index(max(my_y))] , max(my_y))
    plt.text(my_x[my_y.index(max(my_y))] , max(my_y)+5, textstr, bbox=props)

    # Vermits we met dt stapjes werken, komen we niet exact uit op nul. Dit is een benadering van waar het nulpunt ligt.
    nulpunt_x=np.abs(my_x[-2]+(my_x[-1]-my_x[-2])/(my_y[-2]-my_y[-1])*my_y[-2])
    time_of_flight=np.abs(my_t[-2]+(my_t[-1]-my_t[-2])/(my_y[-2]-my_y[-1])*my_y[-2])

    textstr = '{:.0f}m,{:.2f} s'.format(nulpunt_x, time_of_flight)
    # manueel aan te passen indien de labels overlappen.
    if index==3:
        plt.text(max(my_x),-25,textstr,bbox=props)
    else:
        plt.text(max(my_x),5,textstr,bbox=props)
    
    
    # Plot nu ook het traject zonder luchtweerstand in stippellijn
    my_x, my_y, my_t= simuleer_traject(launch_x, launch_y, ss[index], hoek, dt, num_timesteps, dd[index] , CC[index], rr[index], mm[index],False,stop_y)   
    plt.plot(my_x, my_y,color=my_colors[index],linestyle='dotted' )

plt.grid(True)
    
# Bewaar het bestand
dpi=300
myfilename='Pijlenboog_traject_van_een_pijl_met_luchtweerstand_'+LANG+'_'+str(dpi)+'dpi.png'
plt.savefig(myfilename,dpi=300, bbox_inches='tight',transparent=True) 
plt.show()